Polysubstance and tr. completion
Analyze the association between polysubstance at admission and tr. compleiton longitudinally along treatments, accounting for irregular observations.
Setup
Code
rm(list = ls())
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
"E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
"G:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")
load(paste0(folder_path,"an_grant_23_24_3.RData"))Packages
Code
system("fc-cache -f -v")[1] 0Code
if(!require(pacman)){install.packages("pacman");library(pacman)}
pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
if (getRversion() != "4.1.2") { stop("Requires R version 4.1.2. Actual: ", getRversion()) }
if(!require(geepack)){install.packages("geepack");library(geepack)}
if(!require(geeM)){install.packages("geeM");library(geeM)}
if(!require(tidyverse)){install.packages("tidyverse");library(tidyverse)}
if(!require(MASS)){install.packages("MASS");library(MASS)}
if(!require(geeasy)){install.packages("geeasy");library(geeasy)}
if(!require(MuMIn)){remotes::install_version("MuMIn", "1.46.0");library(MuMIn)}
if(!require(tableone)){install.packages("tableone");library(tableone)}
if(!require(knitr)){install.packages("knitr");library(knitr)}
if(!require(emmeans)){install.packages("emmeans");library(emmeans)}
if(!require(biostat3)){install.packages("biostat3");library(biostat3)}
if(!require(rms)){install.packages("rms");library(rms)}
if(!require(meta)){install.packages("meta");library(meta)}
try(if(!require(dmetar)){install.packages("dmetar");library(dmetar)})Error : package 'dmetar' is not availableCode
if(!require(metafor)){install.packages("metafor");library(metafor)}
try(if(!require(extrafont)){install.packages("extrafont");library(extrafont)})
try(if(!require(showtext)){install.packages("showtext");library(showtext)})
try(if(!require(mice)){install.packages("mice");library(mice)})
try(if(!require(htmltools)){install.packages("htmltools");library(htmltools)})
if(!require(chisq.posthoc.test)){devtools::install_github("ebbertd/chisq.posthoc.test")}
if(!require(DescTools)){install.packages("DescTools")}
if(!require(vcd)){install.packages("vcd")}
try(extrafont::font_import(pattern = "Times New Roman", prompt = FALSE))Error in data.frame(fontfile = ttfiles, FontName = "", stringsAsFactors = FALSE) :
arguments imply differing number of rows: 0, 1Code
extrafont::loadfonts(device="win")
# Agregar la ruta a la fuente "Times New Roman"
font_add("Times New Roman", "C:/Windows/Fonts/times.ttf")
# Habilitar showtext para gráficos
showtext::showtext_auto()
#https://gist.github.com/avallecam/56af06f46e5544c3af0f46344df20989
#https://stackoverflow.com/questions/49089476/any-updates-to-model-negative-binomial-distribution-data-with-gee-in-r
#https://www.rdocumentation.org/packages/geepack/versions/1.3.4/topics/geese
#https://stackoverflow.com/questions/13946540/negative-binomial-in-gee?rq=1
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
summary2.geem <- function(object, exponentiate = FALSE, digits = 2, ...) {
if (!is.list(object) || !"beta" %in% names(object)) {
stop("Invalid object: Expected a model object with beta coefficients.")
}
# Initialize coefficients matrix
Coefs <- matrix(NA, nrow = length(object$beta), ncol = 8)
rownames(Coefs) <- object$coefnames
colnames(Coefs) <- c("Term", "Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper", "p.value")
Coefs[, "Estimates"] <- object$beta
# Handling the variance calculations
naive <- is.null(object$var) || is.character(object$var) || any(diag(object$var) < 0)
if (naive) {
warning("Robust variance estimate not available or invalid. Using model-based SE.")
Coefs[, "Robust SE"] <- sqrt(diag(object$naiv.var))
} else {
Coefs[, "Robust SE"] <- sqrt(diag(object$var))
}
Coefs[, "Model SE"] <- sqrt(diag(object$naiv.var)) # Model-based SE
# Calculate Wald statistics
Coefs[, "Wald"] <- Coefs[, "Estimates"] / Coefs[, "Robust SE"]
# 95% CI calculations
z_value <- qnorm((1 + 0.95) / 2)
Coefs[, "CI.lower"] <- Coefs[, "Estimates"] - z_value * Coefs[, "Robust SE"]
Coefs[, "CI.upper"] <- Coefs[, "Estimates"] + z_value * Coefs[, "Robust SE"]
# P-value calculation
Coefs[, "p.value"] <- round(2 * pnorm(-abs(Coefs[, "Wald"])), digits = 4)
# Exponentiate estimates and CIs if specified
if (exponentiate) {
Coefs[, c("Estimates", "CI.lower", "CI.upper")] <- exp(Coefs[, c("Estimates", "CI.lower", "CI.upper")])
}
# Round numerical output
decimal_format <- sprintf("%%1.%df", digits)
Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")] <- round(Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")], digits = digits)
Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")] <- sapply(Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")], function(col) {
as.character(sprintf(decimal_format, as.numeric(col)))
})
# Convert matrix to data frame and keep row names
result_df <- as.data.frame(Coefs)
result_df[, "Term"] <- rownames(Coefs)
row.names(result_df) <- NULL
return(result_df)
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
tidy_geese_model <- function(geese_model, conf.level = 0.95) {
# Extract coefficients and their names
coefficients <- geese_model$beta
coefficient_names <- names(coefficients)
# Extract variance of beta coefficients and compute standard errors
variance_beta <- geese_model$vbeta
std_errors <- sqrt(diag(variance_beta))
# Calculate z-values and p-values
z_values <- coefficients / std_errors
p_values <- 2 * pnorm(-abs(z_values), lower.tail = TRUE)
# Calculate confidence intervals
alpha <- 1 - conf.level
z <- qnorm(1 - alpha / 2)
lower_ci <- coefficients - z * std_errors
upper_ci <- coefficients + z * std_errors
# Exponentiate coefficients and confidence intervals for odds ratios
exp_coefficients <- exp(coefficients)
exp_lower_ci <- exp(lower_ci)
exp_upper_ci <- exp(upper_ci)
# Create a tidy data frame
tidy_data <- tibble::tibble(
Term = coefficient_names,
Estimate = exp_coefficients,
Std.Error = std_errors,
z.value = z_values,
p.value = p_values,
CI.Lower = exp_lower_ci,
CI.Upper = exp_upper_ci
) %>% dplyr::mutate_at(c("Estimate", "Std.Error", "z.value", "CI.Lower","CI.Upper"),~sprintf("%1.2f",.))%>% dplyr::mutate_at(c("p.value"),~sprintf("%1.4f",.))
return(tidy_data)
}
options(knitr.kable.NA = '')
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
knitr::knit_hooks$set(time_it = local({
now <- NULL
function(before, options) {
if (before) {
# record the current time before each chunk
now <<- Sys.time()
} else {
# calculate the time difference after a chunk
res <- ifelse(difftime(Sys.time(), now)>(60^2),difftime(Sys.time(), now)/(60^2),difftime(Sys.time(), now)/(60^1))
# return a character string to show the time
x<-ifelse(difftime(Sys.time(), now)>(60^2),paste("Time for this code chunk to run:", round(res,1), "hours"),paste("Time for this code chunk to run:", round(res,1), "minutes"))
paste('<div class="message">', gsub('##', '\n', x),'</div>', sep = '\n')
}
}
}))
knitr::opts_chunk$set(time_it = TRUE)
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
format_cells <- function(df, rows ,cols, value = c("italics", "bold", "strikethrough")){
# select the correct markup
# one * for italics, two ** for bold
map <- setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough"))
markup <- map[value]
for (r in rows){
for(c in cols){
# Make sure values are not factors
df[[c]] <- as.character( df[[c]])
# Update formatting
df[r, c] <- ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup))
}
}
return(df)
}
#To produce line breaks in messages and warnings
knitr::knit_hooks$set(
error = function(x, options) {
paste('\n\n<div class="alert alert-danger">',
gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),
'</div>', sep = '\n')
},
warning = function(x, options) {
paste('\n\n<div class="alert alert-warning">',
gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),
'</div>', sep = '\n')
},
message = function(x, options) {
paste('<div class="message">',
gsub('##', '\n', x),
'</div>', sep = '\n')
}
)
#_#_#_#_#_#_#_#_#_#_#_#_#_
invisible("Function to format CreateTableOne into a database")
as.data.frame.TableOne <- function(x, ...) {capture.output(print(x,
showAllLevels = TRUE, varLabels = T,...) -> x)
y <- as.data.frame(x)
y$characteristic <- dplyr::na_if(rownames(x), "")
y <- y %>%
fill(characteristic, .direction = "down") %>%
dplyr::select(characteristic, everything())
rownames(y) <- NULL
y}
#_#_#_#_#_#_#_#_#_#_#_#_#_
# Austin, P. C. (2009). The Relative Ability of Different Propensity
# Score Methods to Balance Measured Covariates Between
# Treated and Untreated Subjects in Observational Studies. Medical
# Decision Making. https://doi.org/10.1177/0272989X09341755
smd_bin <- function(x,y){
z <- x*(1-x)
t <- y*(1-y)
k <- sum(z,t)
l <- k/2
return((x-y)/sqrt(l))
}
#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_
counting_process_output<-
function(object=NULL){
broom::tidy(object$m, exponentiate=T, conf.int=T) %>%
dplyr::mutate(across(c("estimate","std.error","robust.se","statistic","conf.low","conf.high"),~sprintf("%1.2f",.))) %>%
dplyr::mutate(p.value= sprintf("%1.4f",p.value))
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
tidy_coxph <- function(model, conf.level = 0.95) {
# Extract the coefficients (a named numeric vector)
coefs <- model$coefficients
# Compute robust standard errors from the robust variance matrix
robust_se <- sqrt(diag(model$var))
# Compute z statistics and two-sided p-values
z <- coefs / robust_se
p_value <- 2 * (1 - pnorm(abs(z)))
# Compute the critical z value for the specified confidence level
z_crit <- qnorm(1 - (1 - conf.level) / 2)
# Exponentiate coefficients to get hazard ratios
hazard_ratio <- exp(coefs)
# Compute the lower and upper confidence intervals (exponentiated)
lower_ci <- exp(coefs - z_crit * robust_se)
upper_ci <- exp(coefs + z_crit * robust_se)
# Build and return a tidy data frame
result <- data.frame(
term = names(coefs),
coef = coefs,
robust_se = robust_se,
hazard_ratio = hazard_ratio,
lower_ci = lower_ci,
upper_ci = upper_ci,
z = z,
p_value = p_value,
stringsAsFactors = FALSE
)
return(result)
}
format_hr_ci <- function(tidy_table, digits = 2) {
# Build a format string dynamically using the desired number of digits.
fmt <- paste0("%.", digits, "f (%.", digits, "f, %.", digits, "f)")
# Create the new column using sprintf to format each row.
tidy_table$hr_ci <- with(tidy_table, sprintf(fmt, hazard_ratio, lower_ci, upper_ci))
return(tidy_table)
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
if(.Platform$OS.type == "windows") withAutoprint({
memory.size()
memory.size(TRUE)
memory.limit()
})> memory.size()
[1] 7546.27
> memory.size(TRUE)
[1] 7546.88
> memory.limit()
[1] 32563Code
memory.limit(size=56000)[1] 56000Code
invisible("At the beginning is in wide format")
invisible("At the beginning is in wide format")
nrow(Base_fiscalia_v15f_grant_23_24_long)
#[1] 109756
length(unique(Base_fiscalia_v15f_grant_23_24_long$hash_key))
#[1] 85048
Base_fiscalia_v15f_grant_23_24_long |>
dplyr::filter(motivodeegreso_mod_imp %in% c("En curso", "Muerte", "Derivación")) |>
(\(df) {
print(message(paste0("Dead, Referrals, ongoing treatments, Entries: ", nrow(df))))
print(message(paste0("Dead, Referrals, ongoing treatments, HASHs: ", nrow(distinct(df,hash_key)))))
})()
Base_fiscalia_v15f_grant_23_24_long |>
dplyr::filter(!motivodeegreso_mod_imp %in% c("En curso", "Muerte", "Derivación")) |>
(\(df) {
print(message(paste0("Dead, Referrals, ongoing treatments, Entries: ", nrow(df))))
print(message(paste0("Dead, Referrals, ongoing treatments, HASHs: ", nrow(distinct(df,hash_key)))))
})()
length(unlist(unique(Base_fiscalia_v15f_grant_23_24_long2$hash_key)))
#[exclude ongoing treatments and external referrals]
length(unique(Base_fiscalia_v15f_grant_23_24$hash_key))
length(unique(Base_fiscalia_v15f_grant_23_24_long2$hash_key))
#72404
nrow(Base_fiscalia_v15f_grant_23_24_long2)
#[1] 90075
#
Base_fiscalia_v15f_grant_23_24_long2 |>
dplyr::filter(!hash_key %in% Base_fiscalia_v15f_grant_23_24_long2_miss_proc_multtr$hash_key
) |>
(\(df) {
print(message(paste0("Only one treatment: ", nrow(df))))
print(message(paste0("Only one treatment: ", nrow(distinct(df,hash_key)))))
})()Code
library(DiagrammeR)
gr<-
grViz("
digraph flowchart {
graph [layout = dot, rankdir = TB]
# General node styling
node [fontname = Times, shape = rectangle, fontsize = 10, style = filled, fillcolor = transparent]
# Main flow nodes
original [label = 'Original C1 Dataset\\n(n = 163,146;\\nPatients = 85,722)', fillcolor = lightgray]
c1_dataset [label = 'C1 Dataset\\n(n = 109,756;\\nPatients = 85,048)']
after_discard [label = 'C1 Dataset\\n(n = 90,075;\\nPatients = 72,404)']
final_dataset [label = 'Final C1 Dataset\\n(n = 30,988;\\nPatients = 13,317)', fillcolor = lightgray]
# Discard nodes (aligned between main flow steps)
discard_referrals [label = '•Discard external referrals or ongoing treatments\\l(n = 19,681; Patients = 18,450)\\l']
discard_duplicates [label = '•Discard duplicated entries\\l•Intermediate events of treatment (continuous referrals)\\l']
discard_single [label = '•Discard admitted to only one treatment(Patients = 59,087)\\l']
# Invisible vertices for middle line
v1 [shape = point, width = 0, style = invis]
v2 [shape = point, width = 0, style = invis]
v3 [shape = point, width = 0, style = invis]
# Main flow edges (vertical line)
original -> v1 [arrowhead = none]
v1 -> c1_dataset
c1_dataset -> v2 [arrowhead = none]
v2 -> after_discard
after_discard -> v3 [arrowhead = none]
v3 -> final_dataset
# Discard connections (from the middle line)
v1 -> discard_referrals
v2 -> discard_duplicates
v3 -> discard_single
# Alignment
{ rank = same; discard_referrals; v1 }
{ rank = same; discard_duplicates; v2 }
{ rank = same; discard_single; v3 }
}
",
width = 800,
height = 600)
grCode
unlink(paste0(folder_path, "_doc/_flowchart_files"), recursive = TRUE)
htmlwidgets::saveWidget(gr, paste0(folder_path, "_doc/_flowchart.html"))
webshot::webshot(paste0(folder_path, "_doc/_flowchart.html"),
paste0(folder_path, "_doc/_flowchart.png"),
vwidth = 300*1.2, vheight = 300, zoom=10, expand=100) # Prueba con diferentes coordenadas top, left, width, and height.
0.a. Descriptives
Code
invisible("2024-05-21: test different")
data_mine_miss_restr_proc2exp<-
data_mine_miss_restr_proc2 %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
#"susinidumrec_otr", "susinidum_coc", "susinidum_pbc", "susinidum_mar",
#janitor::tabyl(n_off_acq)
#(Post-Treatment)
# n_post_off_acq
# n_post_off_vio
# n_post_off_sud
# n_post_off_ot
# n_post_off
invisible("To account for variability by treatment setting, we stratified the analysis by setting: basic ambulatory GP intensive ambulatory GP residential WO intensive ambulatory WO residential")
table(subset(Base_fiscalia_v15f_grant_23_24_long2_miss_proc_multtr, is_first_occurrence==1, select= "tipo_de_plan_2_mod"))
basic ambulatory GP intensive ambulatory GP residential
4360 4998 2178
WO intensive ambulatory WO residential
745 1036 Code
data_mine_miss_restr_proc2exp %>%
group_by(hash_key) %>%
summarise(min(fech_ing_num)) %>% nrow()[1] 13317Code
#13317
subset(data_mine_miss_restr_proc2 , is_first_occurrence==1) %>% nrow()[1] 13317Code
#13317
data_mine_miss_restr_proc2_baseline<-
subset(data_mine_miss_restr_proc2exp , is_first_occurrence==1)
variables_vector <- c("tr_outcome", "policonsumo2", "comp_bpsc_y3_severe",
"less_90d_tr1", "log_dias_treat_imp_sin_na",
"edad_al_ing_1", "ano_nac_corr", "susinidumrec_coc",
"susinidumrec_pbc", "susinidumrec_mar", "susinidumrec_otr",
"psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day",
"cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindumrec_coc",
"susprindumrec_pbc", "susprindumrec_mar", "susprindumrec_otr")#,"lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "lag_tr_outcome")Code
attr(data_mine_miss_restr_proc2_baseline$tr_outcome,"label") <- "Complete status of treatment (binary)"
attr(data_mine_miss_restr_proc2_baseline$tr_outcome,"label") <- "Complete status of treatment (binary)"
attr(data_mine_miss_restr_proc2_baseline$policonsumo2,"label") <- "Polysubstance use"
attr(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe,"label") <- "Biopsychosocial compromise (Severe)"
attr(data_mine_miss_restr_proc2_baseline$less_90d_tr1,"label") <- "Treatment duration (binary) (<90 days)"
attr(data_mine_miss_restr_proc2_baseline$log_dias_treat_imp_sin_na,"label") <- "Treatment duration (log-scaled days)"
attr(data_mine_miss_restr_proc2_baseline$edad_al_ing_1,"label") <- "Age at admission to treatment"
attr(data_mine_miss_restr_proc2_baseline$ano_nac_corr,"label") <- "Birth year"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_coc,"label") <- "Primary substance (initial diagnosis): cocaine hydrochloride"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_pbc,"label") <- "Primary substance (initial diagnosis): cocaine base paste"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_mar,"label") <- "Primary substance (initial diagnosis): marijuana"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_otr,"label") <- "Primary substance (initial diagnosis): other"
attr(data_mine_miss_restr_proc2_baseline$psycom_dum_study,"label") <- "Psychiatric comorbidity (ICD-10): In study"
attr(data_mine_miss_restr_proc2_baseline$psycom_dum_with,"label") <- "Psychiatric comorbidity (ICD-10): Diagnosed"
attr(data_mine_miss_restr_proc2_baseline$freq_cons_dum_5day,"label") <- "Daily frequence of primary substance use at admission"
attr(data_mine_miss_restr_proc2_baseline$cond_oc_dum_2inact,"label") <- "Occupational Status: Inactive"
attr(data_mine_miss_restr_proc2_baseline$cond_oc_dum_3unemp,"label") <- "Occupational Status: Unemployed"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_coc,"label") <- "Primary substance at admission to treatment: cocaine hydrochloride"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_pbc,"label") <- "Primary substance at admission to treatment: cocaine base paste"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_mar,"label") <- "Primary substance at admission to treatment: marijuana"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_otr,"label") <- "Primary substance at admission to treatment: other"
attr(data_mine_miss_restr_proc2_baseline$tipo_de_plan_2_mod,"label") <- "Treatment setting"
invisible("Table for imputed dataset with selected covariates")
library(tableone)
vars <- setdiff(c(variables_vector,"tipo_de_plan_2_mod"), "policonsumo2")
factorVars <- setdiff(variables_vector, c("policonsumo2", "edad_al_ing_1", "log_dias_treat_imp_sin_na", "ano_nac_corr"))
tableOne <- CreateTableOne(vars = vars, data = data_mine_miss_restr_proc2_baseline,
factorVars = factorVars, strata = "policonsumo2",
addOverall = TRUE, # <-- This is the key addition
includeNA = TRUE,
test = TRUE) # Remove statistical tests for speed
print(tableOne, smd = TRUE) # Calculate SMD separately if needed Stratified by policonsumo2
Overall 0
n 13317 2383
tr_outcome = 1 (%) 10448 (78.5) 1833 (76.9)
comp_bpsc_y3_severe = 1 (%) 5496 (41.3) 690 (29.0)
less_90d_tr1 = 1 (%) 3269 (24.5) 567 (23.8)
log_dias_treat_imp_sin_na (mean (SD)) 5.00 (1.07) 5.03 (0.93)
edad_al_ing_1 (mean (SD)) 33.80 (9.46) 38.29 (11.32)
ano_nac_corr (mean (SD)) 1979.02 (9.54) 1975.34 (11.02)
susinidumrec_coc = 1 (%) 627 ( 4.7) 108 ( 4.5)
susinidumrec_pbc = 1 (%) 1123 ( 8.4) 251 (10.5)
susinidumrec_mar = 1 (%) 4348 (32.6) 483 (20.3)
susinidumrec_otr = 1 (%) 295 ( 2.2) 57 ( 2.4)
psycom_dum_study = 1 (%) 2653 (19.9) 420 (17.6)
psycom_dum_with = 1 (%) 5836 (43.8) 986 (41.4)
freq_cons_dum_5day = 1 (%) 6242 (46.9) 1013 (42.5)
cond_oc_dum_2inact = 1 (%) 2410 (18.1) 468 (19.6)
cond_oc_dum_3unemp = 1 (%) 5271 (39.6) 764 (32.1)
susprindumrec_coc = 1 (%) 2370 (17.8) 292 (12.3)
susprindumrec_pbc = 1 (%) 6898 (51.8) 902 (37.9)
susprindumrec_mar = 1 (%) 717 ( 5.4) 64 ( 2.7)
susprindumrec_otr = 1 (%) 188 ( 1.4) 43 ( 1.8)
tipo_de_plan_2_mod (%)
basic ambulatory 4360 (32.7) 1040 (43.6)
GP intensive ambulatory 4998 (37.5) 786 (33.0)
GP residential 2178 (16.4) 272 (11.4)
WO intensive ambulatory 745 ( 5.6) 138 ( 5.8)
WO residential 1036 ( 7.8) 147 ( 6.2)
Stratified by policonsumo2
1 p test SMD
n 10934
tr_outcome = 1 (%) 8615 (78.8) 0.047 0.045
comp_bpsc_y3_severe = 1 (%) 4806 (44.0) <0.001 0.315
less_90d_tr1 = 1 (%) 2702 (24.7) 0.359 0.021
log_dias_treat_imp_sin_na (mean (SD)) 4.99 (1.10) 0.074 0.043
edad_al_ing_1 (mean (SD)) 32.82 (8.70) <0.001 0.542
ano_nac_corr (mean (SD)) 1979.83 (8.99) <0.001 0.446
susinidumrec_coc = 1 (%) 519 ( 4.7) 0.693 0.010
susinidumrec_pbc = 1 (%) 872 ( 8.0) <0.001 0.088
susinidumrec_mar = 1 (%) 3865 (35.3) <0.001 0.341
susinidumrec_otr = 1 (%) 238 ( 2.2) 0.569 0.014
psycom_dum_study = 1 (%) 2233 (20.4) 0.002 0.071
psycom_dum_with = 1 (%) 4850 (44.4) 0.008 0.060
freq_cons_dum_5day = 1 (%) 5229 (47.8) <0.001 0.107
cond_oc_dum_2inact = 1 (%) 1942 (17.8) 0.033 0.048
cond_oc_dum_3unemp = 1 (%) 4507 (41.2) <0.001 0.191
susprindumrec_coc = 1 (%) 2078 (19.0) <0.001 0.187
susprindumrec_pbc = 1 (%) 5996 (54.8) <0.001 0.346
susprindumrec_mar = 1 (%) 653 ( 6.0) <0.001 0.162
susprindumrec_otr = 1 (%) 145 ( 1.3) 0.090 0.039
tipo_de_plan_2_mod (%) <0.001 0.298
basic ambulatory 3320 (30.4)
GP intensive ambulatory 4212 (38.5)
GP residential 1906 (17.4)
WO intensive ambulatory 607 ( 5.6)
WO residential 889 ( 8.1) Code
smd_bin(
prop.table(table(data_mine_miss_restr_proc2_baseline$susprindum_mar[data_mine_miss_restr_proc2_baseline$policonsumo==0]))[[2]],
prop.table(table(data_mine_miss_restr_proc2_baseline$susprindum_mar[data_mine_miss_restr_proc2_baseline$policonsumo==1]))[[2]]
)[1] -0.1630994Code
smd_bin(
prop.table(table(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe [data_mine_miss_restr_proc2_baseline$policonsumo==0]))[[2]],
prop.table(table(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe [data_mine_miss_restr_proc2_baseline$policonsumo==1]))[[2]]
)[1] -0.3140315Code
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
"E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
"C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")
as.data.frame.TableOne(tableOne, smd=T, nonnormal= T)%>%
dplyr::mutate(char2=characteristic) %>%
tidyr::fill(char2) %>%
dplyr::select(char2,everything()) %>%
dplyr::mutate(level=ifelse(is.na(level),"[Missing]",level)) %>%
dplyr::mutate(char2=dplyr::case_when(characteristic=="NA"~NA_character_,TRUE~as.character(characteristic))) %>%
format_cells(1, 1:length(names(.)), "bold") %>%
dplyr::select(-1) %>%
dplyr::mutate_all(~ str_replace_all(., pattern = "\\( ", replacement = "\\(")) %>%
dplyr::mutate_all(~ str_trim(.)) %>%
dplyr::select(characteristic, level, `0`, `1`, Overall, SMD)%>%
{
knitr::kable(.,size=10, format="html",caption= "Summary descriptives, Polysubstance(1) and no Polysubstance use (0)", escape=T)%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
write.table(.,paste0(folder_path,"baseline_psu_desc_subsamp_smd.csv"), dec=",", sep="\t")
}Code
rbind.data.frame(
cbind.data.frame(
type= "PSU at admission and at least one dropout from the first admission",
biostat3::survRate(Surv(cens_time, sum_tr_outcomes>0) ~ policonsumo2,
data= data_mine_miss_restr_proc2 %>%
dplyr::group_by(hash_key) %>%
dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T)) %>%
dplyr::filter(is_first_occurrence==1) %>%
dplyr::ungroup()) %>%
dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000)))
),
cbind.data.frame(
type= "PSU at admission and first dropout",
biostat3::survRate(Surv(cens_time, tr_outcome==1) ~ policonsumo2,
data= data_mine_miss_restr_proc2 %>%
dplyr::group_by(hash_key) %>%
dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T)) %>%
dplyr::filter(is_first_occurrence==1) %>%
dplyr::ungroup()) %>%
dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000)))
),
cbind.data.frame(
type= "At least one treatment reporting PSU and at least one dropout from the first admission",
biostat3::survRate(Surv(cens_time, sum_tr_outcomes>0) ~ policonsumo2,
data= data_mine_miss_restr_proc2 %>%
dplyr::group_by(hash_key) %>%
dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T),
sum_policonsumo2=sum(policonsumo2==1,na.rm=T)) %>%
dplyr::filter(is_first_occurrence==1) %>%
dplyr::ungroup() %>% dplyr::mutate(policonsumo2=ifelse(sum_policonsumo2>0,1,0))) %>%
dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000)))
),
cbind.data.frame(
type= "At least one treatment reporting PSU and first dropout",
biostat3::survRate(Surv(cens_time, tr_outcome==1) ~ policonsumo2,
data= data_mine_miss_restr_proc2 %>%
dplyr::group_by(hash_key) %>%
dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T),
sum_policonsumo2=sum(policonsumo2==1,na.rm=T)) %>%
dplyr::filter(is_first_occurrence==1) %>%
dplyr::ungroup() %>% dplyr::mutate(policonsumo2=ifelse(sum_policonsumo2>0,1,0))) %>%
dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000)))
)
) %>%
dplyr::mutate(`IR (95% CI)`= paste0(rate, " (", lower,", ", upper,")")) %>%
dplyr::select(-any_of(c("rate", "lower", "upper"))) %>%
dplyr::mutate(across(c("tstop", "event"),~ format(as.numeric(.), big.mark=","))) %>%
{
copiar_nombres2(.)
write.table(.,paste0(folder_path,"person_time_irrs.csv"), dec=",", sep="\t")
knitr::kable(.,"html", caption= "Incidence rates (per 1.000 person-months)")%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| type | policonsumo2 | tstop | event | IR (95% CI) | |
|---|---|---|---|---|---|
| policonsumo2=0 | PSU at admission and at least one dropout from the first admission | 0 | 161,852.30 | 2,135 | 13.19 (12.64, 13.76) |
| policonsumo2=1 | PSU at admission and at least one dropout from the first admission | 1 | 872,863.25 | 10,085 | 11.55 (11.33, 11.78) |
| policonsumo2=01 | PSU at admission and first dropout | 0 | 161,852.30 | 1,833 | 11.33 (10.81, 11.86) |
| policonsumo2=11 | PSU at admission and first dropout | 1 | 872,863.25 | 8,615 | 9.87 (9.66, 10.08) |
| policonsumo2=02 | At least one treatment reporting PSU and at least one dropout from the first admission | 0 | 78,926.85 | 1,099 | 13.92 (13.11, 14.77) |
| policonsumo2=12 | At least one treatment reporting PSU and at least one dropout from the first admission | 1 | 955,788.69 | 11,121 | 11.64 (11.42, 11.85) |
| policonsumo2=03 | At least one treatment reporting PSU and first dropout | 0 | 78,926.85 | 936 | 11.86 (11.11, 12.64) |
| policonsumo2=13 | At least one treatment reporting PSU and first dropout | 1 | 955,788.69 | 9,512 | 9.95 (9.75, 10.15) |
0.b. IIWs
Code
rbind.data.frame(
cbind.data.frame(model= "Primary (no corrections, lag=0)", t(matrix(summary(data_mine_miss_restr_proc2$iiw_nocorr_st)))),
cbind.data.frame(model= "Alternative (no corrections, lag=1)", t(matrix(summary(data_mine_miss_restr_proc2$iiw_nocorr_alt_st)))),
cbind.data.frame(model= "Primary (after PH, lag=0)", t(matrix(summary(data_mine_miss_restr_proc2_iiw_after_ph$iiw_after_ph_st)))),
cbind.data.frame(model= "Alternative (after PH, lag=1)", t(matrix(summary(data_mine_miss_restr_proc2_iiw_after_ph_alt$iiw_after_ph_st)))),
cbind.data.frame(model= "Primary (strata, lag=0)", t(matrix(summary(data_mine_miss_proc2_iiw_strata_alt$iiw_strata_st)))),
cbind.data.frame(model= "Alternative (strata, lag=1)", t(matrix(summary(data_mine_miss_proc2_iiw_strata_alt_alt$iiw_strata_st))))
) %>% #Min. First quartile Median Mean Third quartile Max.
dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)) %>%
dplyr::mutate(model= dplyr::case_when(model== "Primary (no corrections, lag=0)"~ "No time-varying coefficients, lagged covariates fixed in 0", model== "Alternative (no corrections, lag=1)"~ "No time-varying coefficients, lagged covariates fixed in 1", model== "Primary (after PH, lag=0)"~ "With time-varying coefficients, lagged covariates fixed in 0", model== "Alternative (after PH, lag=1)"~ "With time-varying coefficients, lagged covariates fixed in 1", model== "Primary (strata, lag=0)"~ "Stratified by follow-up intervals, lagged covariates fixed in 0", model== "Alternative (strata, lag=1)"~ "Stratified by follow-up intervals, lagged covariates fixed in 1")) %>%
{
copiar_nombres2(.)
write.table(.,paste0(folder_path,"iiws_20240516.csv"), dec=",", sep="\t")
knitr::kable(.,size=10, format="html",caption= "Descriptive characterization of inverse intensity weights", escape=T, col.names= c("Visit intensity model", "Min.", "First quartile", "Median", "Mean", "Third quartile", "Max."))%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| Visit intensity model | Min. | First quartile | Median | Mean | Third quartile | Max. |
|---|---|---|---|---|---|---|
| No time-varying coefficients, lagged covariates fixed in 0 | 0.21 | 1.00 | 1.62 | 1.83 | 2.56 | 3.79 |
| No time-varying coefficients, lagged covariates fixed in 1 | 0.13 | 0.42 | 0.68 | 0.70 | 1.00 | 1.00 |
| With time-varying coefficients, lagged covariates fixed in 0 | 0.47 | 0.67 | 1.08 | 1.33 | 1.90 | 3.03 |
| With time-varying coefficients, lagged covariates fixed in 1 | 0.10 | 0.17 | 0.55 | 0.88 | 1.36 | 3.56 |
| Stratified by follow-up intervals, lagged covariates fixed in 0 | 0.29 | 0.60 | 0.86 | 0.99 | 1.18 | 3.02 |
| Stratified by follow-up intervals, lagged covariates fixed in 1 | 0.12 | 0.12 | 0.44 | 0.60 | 0.76 | 3.07 |
0.c. Intensity model
Code
#Cox Counting process
#https://rdrr.io/github/gforge/Greg/src/R/tidy.rms.R
source(paste0(folder_path,"tidy_rms.R"))
term_mapping <- data.frame(
term = c("lag_tr_outcome", "lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "edad_al_ing_1", "ano_nac_corr", "susinidumrec_coc", "susinidumrec_pbc", "susinidumrec_mar", "susinidumrec_otr", "psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day", "cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindumrec_coc", "susprindumrec_pbc", "susprindumrec_mar", "susprindumrec_otr", "lag_tr_outcome_rec", "lag_less_90d_tr1_rec", "lag_comp_bpsc_y3_severe_rec", "susinidum_coc_rec2", "psycom_dum_with_rec2"),
term_label = c("Treatment outcome of the previous treatment", "Previous biopsychosocial compromise (severe)",
"Previous treatment duration (<90 days)", "Previous treatment duration (in logarithmic scaled days)",
"Polysubstance use status of the previous treatment", "Age at admission to treatment", "Birth year",
"Primary substance (initial diagnosis), cocaine", "Primary substance (initial diagnosis), cocaine base paste", "Primary substance (initial diagnosis), marijuana", "Primary substance (initial diagnosis), other", "Psychiatric comorbidity (diagnosis unknown or under study)", "Psychiatric comorbidity (confirmed comorbidity)", "Daily frequence of primary substance use at admission", "Occupational status (inactive)", "Occupational status (unemployed)", "Primary substance at admission to treatment (cocaine hydrochloride)", "Primary substance at admission to treatment (cocaine base paste)", "Primary substance at admission to treatment (marijuana)", "Primary substance at admission to treatment (other)", "Treatment outcome of the previous treatment (recoded for interaction with time)", "Previous biopsychosocial compromise (severe, recoded for interaction with time)", "Previous treatment duration (<90 days, recoded for interaction with time)", "Primary substance (initial diagnosis), cocaine (recoded for interaction with time)", "Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time)")
)
cph_nocorr<-
cph(Surv(lag_time,time,event)~
cluster(id)+
lag_tr_outcome+
lag_comp_bpsc_y3_severe+
lag_less_90d_tr1+
log_lag_dias_treat_imp_sin_na +
lag_policonsumo2 +
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr +
strat(tipo_de_plan_2_mod),
data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(),
x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)
cph_after_ph<-
cph(Surv(lag_time,time,event==1)~
cluster(id)+
lag_tr_outcome_rec +
log_lag_dias_treat_imp_sin_na +
lag_less_90d_tr1_rec+
lag_comp_bpsc_y3_severe_rec +
lag_policonsumo2 +
edad_al_ing_1 +
ano_nac_corr +
susinidum_coc_rec2 +
susinidumrec_otr +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_with_rec2 +
psycom_dum_study +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr +
strat(tipo_de_plan_2_mod),
data=data_mine_miss_restr_proc2exp,
x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)
model_after_time_strat<-
cph(Surv(lag_time,time,event)~
cluster(id)+
lag_tr_outcome+
lag_comp_bpsc_y3_severe+
lag_less_90d_tr1+
log_lag_dias_treat_imp_sin_na +
lag_policonsumo2 +
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr +
strat(tipo_de_plan_2_mod)+
strat(time_interval3),
data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(),
x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)
model_after_time_strat_alt<-
cph(Surv(lag_time,time,event)~
cluster(id)+
lag_tr_outcome+
lag_comp_bpsc_y3_severe+
lag_less_90d_tr1+
log_lag_dias_treat_imp_sin_na +
lag_policonsumo2 +
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_otr +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
strat(tipo_de_plan_2_mod)+
strat(time_interval3_alt),
data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(),
x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)
invisible("Problemático")
# broom::tidy(coxph(Surv(lag_time,time,event==1)~
# cluster(id)+
# lag_tr_outcome_rec +
# log_lag_dias_treat_imp_sin_na +
# lag_less_90d_tr1_rec+
# lag_comp_bpsc_y3_severe_rec +
# lag_policonsumo2 +
# edad_al_ing_1 +
# ano_nac_corr +
# susinidum_coc_rec2 +
# susinidumrec_otr +
# susinidum_pbc +
# susinidum_mar +
# psycom_dum_with_rec2 +
# psycom_dum_study +
# freq_cons_dum_5day +
# cond_oc_dum_2inact +
# cond_oc_dum_3unemp +
# susprindum_oh +
# susprindum_coc +
# susprindum_pbc +
# susprindum_mar+
# strata(tipo_de_plan_2_mod),
# data=data_mine_miss_restr_proc2), exponentiate=T, conf.int=T)
models_hr_table2<-
rbind.data.frame(
cbind.data.frame(model="No correction", tidy.rms(cph_nocorr, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>%
dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="After PH", tidy.rms(cph_after_ph, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>%
dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="Stratifying", tidy.rms(model_after_time_strat, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>%
dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="Stratifying, alt strata", tidy.rms(model_after_time_strat_alt, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")] %>%
dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)))
)
models_hr_table2$term_label <- term_mapping$term_label[match(models_hr_table2$term, term_mapping$term)]
models_hr_table2%>%
dplyr::select(term_label, everything()) %>%
{
#copiar_nombres2(.)
write.table(.,paste0(folder_path,"coxph_intensity_model_prev.csv"), dec=",", sep="\t")
knitr::kable(.,size=10, format="html",caption= "Intensity model", escape=T)%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| term_label | model | term | estimate | conf.low | conf.high |
|---|---|---|---|---|---|
| Treatment outcome of the previous treatment | No correction | lag_tr_outcome | 1.17 | 1.13 | 1.21 |
| Previous biopsychosocial compromise (severe) | No correction | lag_comp_bpsc_y3_severe | 1.06 | 1.03 | 1.10 |
| Previous treatment duration (<90 days) | No correction | lag_less_90d_tr1 | 1.11 | 1.06 | 1.16 |
| Previous treatment duration (in logarithmic scaled days) | No correction | log_lag_dias_treat_imp_sin_na | 0.98 | 0.96 | 1.00 |
| Polysubstance use status of the previous treatment | No correction | lag_policonsumo2 | 0.98 | 0.95 | 1.02 |
| Age at admission to treatment | No correction | edad_al_ing_1 | 1.26 | 1.25 | 1.27 |
| Birth year | No correction | ano_nac_corr | 1.27 | 1.26 | 1.28 |
| Primary substance (initial diagnosis), other | No correction | susinidumrec_otr | 0.96 | 0.87 | 1.05 |
| Primary substance (initial diagnosis), cocaine | No correction | susinidumrec_coc | 1.07 | 0.99 | 1.15 |
| Primary substance (initial diagnosis), cocaine base paste | No correction | susinidumrec_pbc | 1.12 | 1.06 | 1.19 |
| Primary substance (initial diagnosis), marijuana | No correction | susinidumrec_mar | 1.12 | 1.08 | 1.16 |
| Psychiatric comorbidity (diagnosis unknown or under study) | No correction | psycom_dum_study | 1.03 | 0.99 | 1.08 |
| Psychiatric comorbidity (confirmed comorbidity) | No correction | psycom_dum_with | 1.02 | 0.99 | 1.05 |
| Daily frequence of primary substance use at admission | No correction | freq_cons_dum_5day | 1.01 | 0.98 | 1.05 |
| Occupational status (inactive) | No correction | cond_oc_dum_2inact | 1.06 | 1.02 | 1.11 |
| Occupational status (unemployed) | No correction | cond_oc_dum_3unemp | 1.06 | 1.03 | 1.10 |
| Primary substance at admission to treatment (cocaine hydrochloride) | No correction | susprindumrec_coc | 0.99 | 0.95 | 1.04 |
| Primary substance at admission to treatment (cocaine base paste) | No correction | susprindumrec_pbc | 1.01 | 0.97 | 1.05 |
| Primary substance at admission to treatment (marijuana) | No correction | susprindumrec_mar | 0.99 | 0.92 | 1.06 |
| Primary substance at admission to treatment (other) | No correction | susprindumrec_otr | 1.11 | 0.98 | 1.25 |
| Treatment outcome of the previous treatment (recoded for interaction with time) | After PH | lag_tr_outcome_rec | 0.09 | 0.08 | 0.10 |
| Previous treatment duration (in logarithmic scaled days) | After PH | log_lag_dias_treat_imp_sin_na | 0.90 | 0.87 | 0.92 |
| Previous biopsychosocial compromise (severe, recoded for interaction with time) | After PH | lag_less_90d_tr1_rec | 0.62 | 0.49 | 0.77 |
| Previous treatment duration (<90 days, recoded for interaction with time) | After PH | lag_comp_bpsc_y3_severe_rec | 0.19 | 0.16 | 0.22 |
| Polysubstance use status of the previous treatment | After PH | lag_policonsumo2 | 0.99 | 0.95 | 1.03 |
| Age at admission to treatment | After PH | edad_al_ing_1 | 1.20 | 1.19 | 1.21 |
| Birth year | After PH | ano_nac_corr | 1.21 | 1.20 | 1.22 |
| Primary substance (initial diagnosis), cocaine (recoded for interaction with time) | After PH | susinidum_coc_rec2 | 1.04 | 0.94 | 1.15 |
| Primary substance (initial diagnosis), other | After PH | susinidumrec_otr | 1.00 | 0.89 | 1.12 |
| Primary substance (initial diagnosis), cocaine base paste | After PH | susinidumrec_pbc | 1.14 | 1.07 | 1.21 |
| Primary substance (initial diagnosis), marijuana | After PH | susinidumrec_mar | 1.13 | 1.09 | 1.17 |
| Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time) | After PH | psycom_dum_with_rec2 | 1.05 | 1.01 | 1.09 |
| Psychiatric comorbidity (diagnosis unknown or under study) | After PH | psycom_dum_study | 1.17 | 1.12 | 1.23 |
| Daily frequence of primary substance use at admission | After PH | freq_cons_dum_5day | 1.05 | 1.02 | 1.09 |
| Occupational status (inactive) | After PH | cond_oc_dum_2inact | 1.08 | 1.03 | 1.13 |
| Occupational status (unemployed) | After PH | cond_oc_dum_3unemp | 1.10 | 1.06 | 1.14 |
| Primary substance at admission to treatment (cocaine hydrochloride) | After PH | susprindumrec_coc | 0.98 | 0.93 | 1.04 |
| Primary substance at admission to treatment (cocaine base paste) | After PH | susprindumrec_pbc | 1.05 | 1.00 | 1.10 |
| Primary substance at admission to treatment (marijuana) | After PH | susprindumrec_mar | 0.95 | 0.88 | 1.03 |
| Primary substance at admission to treatment (other) | After PH | susprindumrec_otr | 1.00 | 0.87 | 1.14 |
| Treatment outcome of the previous treatment | Stratifying | lag_tr_outcome | 1.15 | 1.11 | 1.19 |
| Previous biopsychosocial compromise (severe) | Stratifying | lag_comp_bpsc_y3_severe | 1.05 | 1.02 | 1.09 |
| Previous treatment duration (<90 days) | Stratifying | lag_less_90d_tr1 | 1.11 | 1.06 | 1.17 |
| Previous treatment duration (in logarithmic scaled days) | Stratifying | log_lag_dias_treat_imp_sin_na | 0.98 | 0.96 | 1.00 |
| Polysubstance use status of the previous treatment | Stratifying | lag_policonsumo2 | 0.98 | 0.94 | 1.02 |
| Age at admission to treatment | Stratifying | edad_al_ing_1 | 1.22 | 1.21 | 1.23 |
| Birth year | Stratifying | ano_nac_corr | 1.22 | 1.21 | 1.23 |
| Primary substance (initial diagnosis), other | Stratifying | susinidumrec_otr | 0.96 | 0.87 | 1.05 |
| Primary substance (initial diagnosis), cocaine | Stratifying | susinidumrec_coc | 1.07 | 0.99 | 1.15 |
| Primary substance (initial diagnosis), cocaine base paste | Stratifying | susinidumrec_pbc | 1.10 | 1.04 | 1.17 |
| Primary substance (initial diagnosis), marijuana | Stratifying | susinidumrec_mar | 1.12 | 1.09 | 1.16 |
| Psychiatric comorbidity (diagnosis unknown or under study) | Stratifying | psycom_dum_study | 1.01 | 0.97 | 1.06 |
| Psychiatric comorbidity (confirmed comorbidity) | Stratifying | psycom_dum_with | 1.01 | 0.98 | 1.05 |
| Daily frequence of primary substance use at admission | Stratifying | freq_cons_dum_5day | 1.01 | 0.98 | 1.04 |
| Occupational status (inactive) | Stratifying | cond_oc_dum_2inact | 1.08 | 1.03 | 1.13 |
| Occupational status (unemployed) | Stratifying | cond_oc_dum_3unemp | 1.09 | 1.05 | 1.13 |
| Primary substance at admission to treatment (cocaine hydrochloride) | Stratifying | susprindumrec_coc | 1.00 | 0.96 | 1.05 |
| Primary substance at admission to treatment (cocaine base paste) | Stratifying | susprindumrec_pbc | 1.02 | 0.98 | 1.06 |
| Primary substance at admission to treatment (marijuana) | Stratifying | susprindumrec_mar | 1.00 | 0.93 | 1.08 |
| Primary substance at admission to treatment (other) | Stratifying | susprindumrec_otr | 1.08 | 0.95 | 1.24 |
| Treatment outcome of the previous treatment | Stratifying, alt strata | lag_tr_outcome | 1.05 | 1.01 | 1.09 |
| Previous biopsychosocial compromise (severe) | Stratifying, alt strata | lag_comp_bpsc_y3_severe | 1.02 | 0.99 | 1.06 |
| Previous treatment duration (<90 days) | Stratifying, alt strata | lag_less_90d_tr1 | 1.04 | 0.99 | 1.09 |
| Previous treatment duration (in logarithmic scaled days) | Stratifying, alt strata | log_lag_dias_treat_imp_sin_na | 0.97 | 0.96 | 0.99 |
| Polysubstance use status of the previous treatment | Stratifying, alt strata | lag_policonsumo2 | 0.99 | 0.95 | 1.03 |
| Age at admission to treatment | Stratifying, alt strata | edad_al_ing_1 | 1.05 | 1.05 | 1.06 |
| Birth year | Stratifying, alt strata | ano_nac_corr | 1.05 | 1.05 | 1.06 |
| Primary substance (initial diagnosis), other | Stratifying, alt strata | susinidumrec_otr | 1.03 | 0.93 | 1.14 |
| Primary substance (initial diagnosis), cocaine | Stratifying, alt strata | susinidumrec_coc | 1.10 | 1.02 | 1.18 |
| Primary substance (initial diagnosis), cocaine base paste | Stratifying, alt strata | susinidumrec_pbc | 1.02 | 0.97 | 1.08 |
| Primary substance (initial diagnosis), marijuana | Stratifying, alt strata | susinidumrec_mar | 1.02 | 0.99 | 1.06 |
| Psychiatric comorbidity (diagnosis unknown or under study) | Stratifying, alt strata | psycom_dum_study | 1.00 | 0.96 | 1.05 |
| Psychiatric comorbidity (confirmed comorbidity) | Stratifying, alt strata | psycom_dum_with | 0.98 | 0.95 | 1.01 |
| Daily frequence of primary substance use at admission | Stratifying, alt strata | freq_cons_dum_5day | 1.02 | 0.99 | 1.05 |
| Occupational status (inactive) | Stratifying, alt strata | cond_oc_dum_2inact | 0.99 | 0.94 | 1.03 |
| Occupational status (unemployed) | Stratifying, alt strata | cond_oc_dum_3unemp | 1.02 | 0.99 | 1.06 |
| Primary substance at admission to treatment (other) | Stratifying, alt strata | susprindumrec_otr | 0.92 | 0.80 | 1.06 |
| Primary substance at admission to treatment (cocaine hydrochloride) | Stratifying, alt strata | susprindumrec_coc | 0.98 | 0.93 | 1.03 |
| Primary substance at admission to treatment (cocaine base paste) | Stratifying, alt strata | susprindumrec_pbc | 0.98 | 0.94 | 1.03 |
| Primary substance at admission to treatment (marijuana) | Stratifying, alt strata | susprindumrec_mar | 0.99 | 0.92 | 1.06 |
0.c.Intensity model (with model that originated IWWs)
Taking other substances as dummy just as in irrelong-21-iiw-weights-nocorr (step 3).
Code
term_mapping2 <- data.frame(
term = c("lag_tr_outcome", "lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "edad_al_ing_1", "ano_nac_corr", "susinidum_coc", "susinidum_pbc", "susinidum_mar", "susinidum_oh", "psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day", "cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh", "lag_tr_outcome_rec", "lag_less_90d_tr1_rec", "lag_comp_bpsc_y3_severe_rec", "susinidum_coc_rec2", "psycom_dum_with_rec2"),
term_label = c("Treatment outcome of the previous treatment", "Previous biopsychosocial compromise (severe)",
"Previous treatment duration (<90 days)", "Previous treatment duration (in logarithmic scaled days)",
"Polysubstance use status of the previous treatment", "Age at admission to treatment", "Birth year",
"Primary substance (initial diagnosis), cocaine", "Primary substance (initial diagnosis), cocaine base paste", "Primary substance (initial diagnosis), marijuana", "Primary substance (initial diagnosis), alcohol", "Psychiatric comorbidity (diagnosis unknown or under study)", "Psychiatric comorbidity (confirmed comorbidity)", "Daily frequence of primary substance use at admission", "Occupational status (inactive)", "Occupational status (unemployed)", "Primary substance at admission to treatment (cocaine hydrochloride)", "Primary substance at admission to treatment (cocaine base paste)", "Primary substance at admission to treatment (marijuana)", "Primary substance at admission to treatment (alcohol)", "Treatment outcome of the previous treatment (recoded for interaction with time)", "Previous biopsychosocial compromise (severe, recoded for interaction with time)", "Previous treatment duration (<90 days, recoded for interaction with time)", "Primary substance (initial diagnosis), cocaine (recoded for interaction with time)", "Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time)")
)
#susinidumrec_oh susprindumrec_oh
cph_nocorr_oh<-
cph(Surv(lag_time,time,event)~
cluster(id)+
lag_tr_outcome+ #tr_outcome.lag
lag_comp_bpsc_y3_severe+ #comp_bpsc_y3_severe.lag
lag_less_90d_tr1+ #less_90d_tr1.lag
log_lag_dias_treat_imp_sin_na + #log_dias_treat_imp_sin_na.lag
lag_policonsumo2 + #policonsumo2.lag
edad_al_ing_1 +
ano_nac_corr +
susinidum_oh+
#susinidumrec_otr +
susinidum_coc +
susinidum_pbc +
susinidum_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindum_oh +
susprindum_coc +
susprindum_pbc +
susprindum_mar +
strat(tipo_de_plan_2_mod),
data= data_mine_miss_restr_proc2 %>% data.table::as.data.table() %>% data.frame(),
x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)
models_hr_table2_oh<-
rbind.data.frame(
cbind.data.frame(model="No correction", tidy.rms(cph_nocorr_oh, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>% dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)))
)
models_hr_table2_oh$term_label <- term_mapping2$term_label[match(models_hr_table2_oh$term, term_mapping2$term)]
models_hr_table2_oh%>%
dplyr::select(term_label, everything()) %>%
{
#copiar_nombres2(.)
write.table(.,paste0(folder_path,"coxph_intensity_model_prev_250216.csv"), dec=",", sep="\t")
knitr::kable(.,size=10, format="html",caption= "Intensity model (no correction, original)", escape=T)%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| term_label | model | term | estimate | conf.low | conf.high |
|---|---|---|---|---|---|
| Treatment outcome of the previous treatment | No correction | lag_tr_outcome | 1.17 | 1.13 | 1.21 |
| Previous biopsychosocial compromise (severe) | No correction | lag_comp_bpsc_y3_severe | 1.06 | 1.03 | 1.10 |
| Previous treatment duration (<90 days) | No correction | lag_less_90d_tr1 | 1.11 | 1.06 | 1.16 |
| Previous treatment duration (in logarithmic scaled days) | No correction | log_lag_dias_treat_imp_sin_na | 0.98 | 0.96 | 1.00 |
| Polysubstance use status of the previous treatment | No correction | lag_policonsumo2 | 0.98 | 0.95 | 1.02 |
| Age at admission to treatment | No correction | edad_al_ing_1 | 1.26 | 1.25 | 1.27 |
| Birth year | No correction | ano_nac_corr | 1.27 | 1.26 | 1.28 |
| Primary substance (initial diagnosis), alcohol | No correction | susinidum_oh | 1.05 | 0.95 | 1.14 |
| Primary substance (initial diagnosis), cocaine | No correction | susinidum_coc | 1.12 | 1.00 | 1.25 |
| Primary substance (initial diagnosis), cocaine base paste | No correction | susinidum_pbc | 1.17 | 1.06 | 1.30 |
| Primary substance (initial diagnosis), marijuana | No correction | susinidum_mar | 1.17 | 1.07 | 1.29 |
| Psychiatric comorbidity (diagnosis unknown or under study) | No correction | psycom_dum_study | 1.03 | 0.99 | 1.08 |
| Psychiatric comorbidity (confirmed comorbidity) | No correction | psycom_dum_with | 1.02 | 0.99 | 1.05 |
| Daily frequence of primary substance use at admission | No correction | freq_cons_dum_5day | 1.01 | 0.98 | 1.05 |
| Occupational status (inactive) | No correction | cond_oc_dum_2inact | 1.06 | 1.02 | 1.11 |
| Occupational status (unemployed) | No correction | cond_oc_dum_3unemp | 1.06 | 1.03 | 1.10 |
| Primary substance at admission to treatment (alcohol) | No correction | susprindum_oh | 0.90 | 0.80 | 1.02 |
| Primary substance at admission to treatment (cocaine hydrochloride) | No correction | susprindum_coc | 0.90 | 0.79 | 1.01 |
| Primary substance at admission to treatment (cocaine base paste) | No correction | susprindum_pbc | 0.91 | 0.81 | 1.03 |
| Primary substance at admission to treatment (marijuana) | No correction | susprindum_mar | 0.89 | 0.78 | 1.02 |
0.d. Cox intensity model behind weights, under two scenarios (lag1 and lag0)
Code
update_terms <- function(tidy_table, term_mapping2) {
# Function to convert tidy_table term to mapping key.
convert_term <- function(term) {
# Check if the term ends with ".lag"
if(grepl("\\.lag$", term)) {
base <- sub("\\.lag$", "", term)
# Special case: if the base is "log_dias_treat_imp_sin_na", insert "lag" after "log_"
if(base == "log_dias_treat_imp_sin_na") {
return("log_lag_dias_treat_imp_sin_na")
} else {
return(paste0("lag_", base))
}
} else {
return(term)
}
}
# Create a new column with the mapping key.
tidy_table$term_key <- sapply(tidy_table$term, convert_term)
# Merge the tidy_table with the term_mapping on the mapping key.
# term_mapping has columns "term" and "term_label".
tidy_table <- merge(tidy_table, term_mapping2, by.x = "term_key", by.y = "term", all.x = TRUE)
# Replace the original term with the label if available.
tidy_table$term <- ifelse(!is.na(tidy_table$term_label), tidy_table$term_label, tidy_table$term)
# Remove temporary columns used for the merge.
tidy_table$term_key <- NULL
tidy_table$term_label <- NULL
return(tidy_table)
}
tidy_table <- tidy_coxph(iiw_nocorr$m)
tidy_table_alt <- tidy_coxph(iiw_nocorr_alt$m)
term_vector <- c(
"Treatment outcome of the previous treatment",
"Previous treatment duration (in logarithmic scaled days)",
"Previous treatment duration (<90 days)",
"Previous biopsychosocial compromise (severe)",
"Polysubstance use status of the previous treatment",
"Age at admission to treatment",
"Birth year",
"Primary substance (initial diagnosis), cocaine",
"Primary substance (initial diagnosis), alcohol",
"Primary substance (initial diagnosis), cocaine base paste",
"Primary substance (initial diagnosis), marijuana",
"Psychiatric comorbidity (confirmed comorbidity)",
"Psychiatric comorbidity (diagnosis unknown or under study)",
"Daily frequence of primary substance use at admission",
"Occupational status (inactive)",
"Occupational status (unemployed)",
"Primary substance at admission to treatment (alcohol)",
"Primary substance at admission to treatment (cocaine hydrochloride)",
"Primary substance at admission to treatment (cocaine base paste)",
"Primary substance at admission to treatment (marijuana)"
)
rbind.data.frame(cbind.data.frame(lag="Lag 0", update_terms(format_hr_ci(tidy_table, digits = 2), term_mapping2)),cbind.data.frame(lag="Lag 1", update_terms(format_hr_ci(tidy_table_alt, digits = 2), term_mapping2)))%>%
dplyr::mutate(term= factor(term, levels=term_vector))%>%
dplyr::arrange(lag, term)%>% #rio::export("clipboard")
{
#copiar_nombres2(.)
write.table(., paste0(folder_path,"coxph_intensity_model_250216.csv"), dec=",", sep="\t")
knitr::kable(tidyr::pivot_wider(dplyr::select(.,"term","lag", "hr_ci"), names_from="lag", values_from="hr_ci"), size=10, format="html",caption= "Specifications of the treatment (visit) process, in Hazard Ratios (HR) and under different lagged variables scenarios (no correction, original)", escape=T)%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| term | Lag 0 | Lag 1 |
|---|---|---|
| Treatment outcome of the previous treatment | 0.65 (0.63, 0.67) | 1.30 (1.26, 1.34) |
| Previous treatment duration (in logarithmic scaled days) | 0.85 (0.84, 0.86) | 1.48 (1.44, 1.53) |
| Previous treatment duration (<90 days) | 0.69 (0.67, 0.72) | 2.63 (2.52, 2.74) |
| Previous biopsychosocial compromise (severe) | 0.88 (0.85, 0.90) | 1.54 (1.50, 1.57) |
| Polysubstance use status of the previous treatment | 0.62 (0.61, 0.64) | 1.26 (1.22, 1.30) |
| Age at admission to treatment | 1.05 (1.05, 1.06) | 1.06 (1.06, 1.06) |
| Birth year | 1.06 (1.06, 1.06) | 1.06 (1.06, 1.06) |
| Primary substance (initial diagnosis), cocaine | 1.01 (0.95, 1.08) | 1.05 (0.98, 1.12) |
| Primary substance (initial diagnosis), alcohol | 1.01 (0.96, 1.07) | 1.04 (0.99, 1.10) |
| Primary substance (initial diagnosis), cocaine base paste | 0.98 (0.92, 1.04) | 1.05 (0.99, 1.11) |
| Primary substance (initial diagnosis), marijuana | 1.04 (0.99, 1.10) | 1.03 (0.98, 1.09) |
| Psychiatric comorbidity (confirmed comorbidity) | 1.00 (0.99, 1.02) | 0.95 (0.94, 0.97) |
| Psychiatric comorbidity (diagnosis unknown or under study) | 1.00 (0.98, 1.03) | 0.74 (0.72, 0.76) |
| Daily frequence of primary substance use at admission | 1.02 (1.01, 1.04) | 0.95 (0.94, 0.97) |
| Occupational status (inactive) | 0.99 (0.96, 1.01) | 0.96 (0.94, 0.99) |
| Occupational status (unemployed) | 1.03 (1.01, 1.05) | 0.98 (0.96, 1.00) |
| Primary substance at admission to treatment (alcohol) | 0.99 (0.93, 1.06) | 1.03 (0.96, 1.09) |
| Primary substance at admission to treatment (cocaine hydrochloride) | 1.11 (1.04, 1.19) | 1.03 (0.97, 1.10) |
| Primary substance at admission to treatment (cocaine base paste) | 1.14 (1.07, 1.22) | 1.01 (0.95, 1.08) |
| Primary substance at admission to treatment (marijuana) | 1.04 (0.96, 1.12) | 1.05 (0.98, 1.12) |
GEE
1. No correction
Code
plan_names <- attr(table(data_mine_miss_restr_proc2$tipo_de_plan_2_mod),"dimnames")[[1]]1.1. No correciton, no weight
Initialize a list to store models
Code
nocorr_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, #origen_ingreso_mod fis_comorbidity_icd_10 dg_trs_cons_sus_or
id = id,
data = current_data,
family = poisson(),
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}1.1.1 No correction, no weight, GEEM Poisson
Initialize a list to store models
Code
geem_nocorr_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar+
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.11.1.2 No correction, no weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(0.2, 16, by = 0.1)Initialize a list to store models and their QICs
Code
geem2_nocorr_nowgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar+
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 16.0: 19838.564776
Best QIC for plan GP intensive ambulatory with x = 16.0: 22686.829620
Best QIC for plan GP residential with x = 16.0: 9979.192804
Best QIC for plan WO intensive ambulatory with x = 16.0: 3487.845045
Best QIC for plan WO residential with x = 16.0: 4909.0672521.2. No correction, weight
Initialize a list to store models
Code
nocorr_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_nocorr_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}1.2.1 No correction, no weight, GEEM Poisson
Initialize a list to store models
Code
geem_nocorr_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_nocorr_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19465.1
QIC for plan GP intensive ambulatory = 22275.9
QIC for plan GP residential = 9833.7
QIC for plan WO intensive ambulatory = 3434.7
QIC for plan WO residential = 4839.11.2.2 No correction, no weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_nocorr_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_nocorr_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19465.056444
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22275.954182
Best QIC for plan GP residential with x = 1588300.0: 9833.688227
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3434.680800
Best QIC for plan WO residential with x = 1588300.0: 4839.094573Check the result for a particular plan
Code
print(geem2_nocorr_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_nocorr_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-9.22 0.02195 -0.03892 -0.002341 0.004535
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.03033 -0.004384 0.09993 0.02421
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1066 0.03428 0.0419 -0.05004
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.005096 0.01384 0.04043 0.002462
susprindumrec_otr
-0.1842
Scale Parameter: 0.4019
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19465.06 1.3. No correction, alt. weight
Initialize a list to store models
Code
nocorr_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_nocorr_alt_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}1.3.1 No correction, alt weight, GEEM Poisson
Initialize a list to store models
Code
geem_nocorr_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_nocorr_alt_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19458.8
QIC for plan GP intensive ambulatory = 22271.9
QIC for plan GP residential = 9827.8
QIC for plan WO intensive ambulatory = 3431.3
QIC for plan WO residential = 4835.61.3.2 No correction, alt weight, GEEM Negative binomial
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_nocorr_alt_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_nocorr_alt_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19458.771969
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22271.862242
Best QIC for plan GP residential with x = 1588300.0: 9827.753103
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3431.276563
Best QIC for plan WO residential with x = 1588300.0: 4835.575200Check the result for a particular plan
Code
print(geem2_nocorr_alt_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_nocorr_alt_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-19.1 0.02437 -0.03728 0.003479 0.009432
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.07601 -0.005822 0.0957 0.02467
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1475 0.04159 0.03779 -0.04085
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.01445 0.01352 0.04788 0.007057
susprindumrec_otr
-0.1853
Scale Parameter: 0.1465
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19458.77 2. After PH
2.1. After PH, no weight
Initialize a list to store models
Code
afterph_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
data_mine_miss_restr_proc2exp_iiw_after_ph <-
data_mine_miss_restr_proc2_iiw_after_ph %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
#weight= current_data$iiw_nocorr_alt_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}2.1.1 After PH, no weight, GEEM Poisson
Initialize a list to store models
Code
geem_afterph_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.12.1.2 After PH, no weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_afterph_nowgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19453.628013
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22261.947187
Best QIC for plan GP residential with x = 1588300.0: 9818.145865
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3424.239894
Best QIC for plan WO residential with x = 1588300.0: 4826.145839Check the result for a particular plan
Code
print(geem2_nocorr_nowgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence")
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-12.32 0.02507 -0.03872 -0.000283 0.006061
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.05421 -0.003447 0.09881 0.02315
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1383 0.04429 0.04057 -0.04427
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.0111 0.01587 0.04783 0.007784
susprindumrec_otr
-0.1904
Scale Parameter: 0.2033
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 16
$best_qic
QIC
19838.56 As expected, QIC is equal, because the covariates and observations were the same
2.2. After PH, weight
Initialize a list to store models
Code
afterph_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_after_ph_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}2.2.1 After PH, weight, GEEM Poisson
Initialize a list to store models
Code
geem_afterph_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindum_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_after_ph_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19479.4
QIC for plan GP intensive ambulatory = 22284.4
QIC for plan GP residential = 9830.3
QIC for plan WO intensive ambulatory = 3445.9
QIC for plan WO residential = 4839.72.2.2 After PH, weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_afterph_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_after_ph_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19479.434224
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22284.436679
Best QIC for plan GP residential with x = 1588300.0: 9830.254882
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3445.890741
Best QIC for plan WO residential with x = 1588300.0: 4839.713265Check the result for a particular plan
Code
print(geem2_afterph_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_after_ph_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-5.745 0.02667 -0.04339 -0.004904 0.00281
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.0327 -0.01607 0.1156 0.02389
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1134 0.0387 0.04731 -0.05169
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.002699 0.01144 0.03865 0.003677
susprindumrec_otr
-0.1873
Scale Parameter: 0.2933
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19479.43 2.3. After PH, alt. weight
Initialize a list to store models
Code
afterph_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
data_mine_miss_restr_proc2exp_iiw_after_ph_alt <-
data_mine_miss_restr_proc2_iiw_after_ph_alt %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_after_ph_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}2.3.1 After PH, alt weight, GEEM Poisson
Initialize a list to store models
Code
geem_afterph_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_after_ph_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19519.4
QIC for plan GP intensive ambulatory = 22324.6
QIC for plan GP residential = 9868.4
QIC for plan WO intensive ambulatory = 3468.6
QIC for plan WO residential = 4861.02.3.2 After PH, alt weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_afterph_alt_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
weight = current_data$iiw_after_ph_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19519.425246
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22324.624851
Best QIC for plan GP residential with x = 1588300.0: 9868.351832
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3468.579298
Best QIC for plan WO residential with x = 1588300.0: 4861.024120Check the result for a particular plan
Code
print(geem2_afterph_alt_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_after_ph_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
0.7589 0.02652 -0.03843 -0.007929 -0.0004416
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.03431 -0.04666 0.1411 0.03227
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.09837 0.04042 0.0413 -0.06383
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.004229 0.03314 0.0559 0.03698
susprindumrec_otr
-0.1908
Scale Parameter: 0.2291
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19519.43 3. Stratifying
3.1. Stratifying, no weight
Initialize a list to store models
Code
strata_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
data_mine_miss_proc2exp_iiw_strata <-
data_mine_miss_proc2_iiw_strata %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
#weight= current_data$iiw_nocorr_alt_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}3.1.1 Stratifying, no weight, GEEM Poisson
Initialize a list to store models
Code
geem_strata_nowgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.13.1.2 Stratifying, no weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_strata_nowgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, # Include relevant predictors
id = id,
data = current_data,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19453.628013
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22261.947187
Best QIC for plan GP residential with x = 1588300.0: 9818.145865
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3424.239894
Best QIC for plan WO residential with x = 1588300.0: 4826.145839Check the result for a particular plan
Code
print(geem2_nocorr_nowgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence")
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-12.32 0.02507 -0.03872 -0.000283 0.006061
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.05421 -0.003447 0.09881 0.02315
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1383 0.04429 0.04057 -0.04427
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.0111 0.01587 0.04783 0.007784
susprindumrec_otr
-0.1904
Scale Parameter: 0.2033
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 16
$best_qic
QIC
19838.56 3.2. Stratifying, weight
Initialize a list to store models
Code
strata_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_strata_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
strata_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}3.2.1 Stratifying, weight, GEEM Poisson
Initialize a list to store models
Code
geem_strata_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight = current_data$iiw_strata_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_strata_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19476.8
QIC for plan GP intensive ambulatory = 22283.0
QIC for plan GP residential = 9839.4
QIC for plan WO intensive ambulatory = 3453.0
QIC for plan WO residential = 4851.13.2.2 Stratifying, weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_strata_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight= current_data$iiw_strata_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_strata_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19476.842070
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22283.006125
Best QIC for plan GP residential with x = 1588300.0: 9839.424427
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3452.975453
Best QIC for plan WO residential with x = 1588300.0: 4851.130363Check the result for a particular plan
Code
print(geem2_strata_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_strata_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-15.18 0.01948 -0.05564 0.0002688 0.0075
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.05049 0.005513 0.1008 0.02269
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1429 0.05033 0.05562 -0.05253
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.01156 -0.001358 0.04337 0.04282
susprindumrec_otr
-0.1754
Scale Parameter: 0.3845
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19476.84 3.3. Stratifying, alt. weight
Initialize a list to store models
Code
strata_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
data_mine_miss_proc2exp_iiw_strata_alt <-
data_mine_miss_proc2_iiw_strata_alt %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_strata_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}3.3.1 Stratifying, alt weight, GEEM Poisson
Initialize a list to store models
Code
geem_strata_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight = current_data$iiw_strata_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19467.2
QIC for plan GP intensive ambulatory = 22269.4
QIC for plan GP residential = 9837.5
QIC for plan WO intensive ambulatory = 3437.3
QIC for plan WO residential = 4841.8Code
#Alternative, better than mainCode
print(geem_strata_alt_wgt_list[[1]])geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = poisson("log"), corstr = "independence", weights = current_data$iiw_strata_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-9.62 0.03745 -0.0414 -0.001121 0.004702
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.09628 -0.01782 0.1065 0.02583
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1531 0.05717 0.03588 -0.03803
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.01523 0.02253 0.04334 0.01237
susprindumrec_otr
-0.2359
Scale Parameter: 0.1917
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993 3.3.2 No correction, alt weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_strata_alt_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight = current_data$iiw_strata_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19467.162599
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22269.414488
Best QIC for plan GP residential with x = 1588300.0: 9837.476338
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3437.325576
Best QIC for plan WO residential with x = 1588300.0: 4841.777059Check the result for a particular plan
Code
print(geem2_strata_alt_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_strata_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-9.62 0.03745 -0.0414 -0.001121 0.004702
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.09628 -0.01782 0.1065 0.02583
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1531 0.05717 0.03588 -0.03803
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.01523 0.02253 0.04334 0.01237
susprindumrec_otr
-0.2359
Scale Parameter: 0.1917
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19467.16 3.4. Stratifying, 2nd alt. weight
Initialize a list to store models
Code
strata_alt_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
data_mine_miss_proc2exp_iiw_strata_alt_alt <-
data_mine_miss_proc2_iiw_strata_alt_alt %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
)
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geese(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
family = poisson(),
weight = current_data$iiw_strata_st,
corstr = "independence",
jack = T)
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}3.4.1 Stratifying, alt weight, GEEM Poisson
Initialize a list to store models
Code
geem_strata_alt_alt_wgt_list <- list()Loop over each unique plan name to fit a model
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset using geem
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight = current_data$iiw_strata_st,
family = poisson("log"),
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
if(model$converged==FALSE){stop("model did not coverge")}
cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
geem_strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}QIC for plan basic ambulatory = 19505.4
QIC for plan GP intensive ambulatory = 22280.0
QIC for plan GP residential = 9885.3
QIC for plan WO intensive ambulatory = 3466.2
QIC for plan WO residential = 4874.2Code
invisible("Alternative, worst model of strata")3.4.2 Stratifying, alt weight, GEEM Negative binomial
Define a sequence of x values
Code
x_values <- seq(158e4, 159e4, by = 100)Initialize a list to store models and their QICs
Code
geem2_strata_alt_alt_wgt_list <- list()Loop over each unique plan name to fit models
Code
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
# Initialize variable to store the best QIC and associated model
best_qic <- Inf
best_model <- NULL
best_x <- NA
# Loop through each x value to find the one with the lowest QIC
for (x in x_values) {
# Fit the GEE model for the current subset using geem with negative binomial
model <- geem(tr_outcome ~ policonsumo2 +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id = id,
data = current_data,
weight = current_data$iiw_strata_st,
family = MASS::negative.binomial(x),
corstr = "independence")
if (model$converged) {
# Calculate QIC for the model
model_qic <- MuMIn::QIC(model)
# Check if this model has a better QIC
if (model_qic < best_qic) {
best_qic <- model_qic
best_model <- model
best_x <- x
}
}
}
# Store the best model and its details in the list
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
geem2_strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
# Optional: Print out some information on progress
cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}Best QIC for plan basic ambulatory with x = 1588300.0: 19505.358265
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22279.986625
Best QIC for plan GP residential with x = 1588300.0: 9885.267664
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3466.220923
Best QIC for plan WO residential with x = 1588300.0: 4874.236495Check the result for a particular plan
Code
print(geem2_strata_alt_alt_wgt_list[[1]])$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe +
edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc +
susinidumrec_pbc + susinidumrec_mar + psycom_dum_study +
psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact +
cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc +
susprindumrec_mar + susprindumrec_otr, id = id, data = current_data,
family = MASS::negative.binomial(x), corstr = "independence",
weights = current_data$iiw_strata_st)
Coefficients:
(Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
-12.83 0.02767 -0.05236 0.001664 0.00628
susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
0.1033 -0.0514 0.1021 0.01249
psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
0.1768 0.06588 0.02143 -0.03881
cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
0.02918 0.01468 0.04859 0.005722
susprindumrec_otr
-0.2477
Scale Parameter: 0.1097
Correlation Model: independence
Estimated Correlation Parameter: 0
Number of clusters: 4360 Maximum cluster size: 10
Number of observations with nonzero weight: 9993
$best_x
[1] 1588300
$best_qic
QIC
19505.36 Model output
Code
types_of_model<-
c("No time-varying coefficients, no weight",
"No time-varying coefficients, lag=0",
"No time-varying coefficients, lag=1",
"With time-varying coefficients, no weight",
"With time-varying coefficients, lag=0",
"With time-varying coefficients, lag=1",
"Stratified by follow-up intervals, no weight",
"Stratified by follow-up intervals, lag=0",
"Stratified by follow-up intervals, lag=1",
"Stratified by follow-up intervals, 2nd, lag=0"
)
types_of_model2<-
c("No time-varying coefficients, no weight",
"No time-varying coefficients, no weight, NB",
"No time-varying coefficients, lag=0",
"No time-varying coefficients, lag=0, NB",
"No time-varying coefficients, lag=1",
"No time-varying coefficients, lag=1, NB",
"With time-varying coefficients, no weight",
"With time-varying coefficients, no weight, NB",
"With time-varying coefficients, lag=0",
"With time-varying coefficients, lag=0, NB",
"With time-varying coefficients, lag=1",
"With time-varying coefficients, lag=1, NB",
"Stratified by follow-up intervals, no weight",
"Stratified by follow-up intervals, no weight, NB",
"Stratified by follow-up intervals, lag=0",
"Stratified by follow-up intervals, lag=0, NB",
"Stratified by follow-up intervals, lag=1",
"Stratified by follow-up intervals, lag=1, NB",
"Stratified by follow-up intervals, 2nd, lag=1",
"Stratified by follow-up intervals, 2nd, lag=1, NB"
)
plan_names[1]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[1]]), Term=="policonsumo2")))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geese_", gsub(" ","_",plan_names[1])), ., envir = .GlobalEnv)
print(.)
}
cbind.data.frame(model= types_of_model2,
QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[1]])),
as.numeric(geem2_nocorr_nowgt_list[[1]]$best_qic),
as.numeric(QIC(geem_nocorr_wgt_list[[1]])),
as.numeric(geem2_nocorr_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_nocorr_alt_wgt_list[[1]])),
as.numeric(geem2_nocorr_alt_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_afterph_nowgt_list[[1]])),
as.numeric(geem2_afterph_nowgt_list[[1]]$best_qic),
as.numeric(QIC(geem_afterph_wgt_list[[1]])),
as.numeric(geem2_afterph_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_afterph_alt_wgt_list[[1]])),
as.numeric(geem2_afterph_alt_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_strata_nowgt_list[[1]])),
as.numeric(geem2_strata_nowgt_list[[1]]$best_qic),
as.numeric(QIC(geem_strata_wgt_list[[1]])),
as.numeric(geem2_strata_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_strata_alt_wgt_list[[1]])),
as.numeric(geem2_strata_alt_wgt_list[[1]]$best_qic),
as.numeric(QIC(geem_strata_alt_alt_wgt_list[[1]])),
as.numeric(geem2_strata_alt_alt_wgt_list[[1]]$best_qic)
),
rbind.data.frame(
dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geem_", gsub(" ","_",plan_names[1])), ., envir = .GlobalEnv)
print(.)
}
i=2
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2")))%>%
{
copiar_nombres2(.)
print(.)
}
#[1] "GP intensive ambulatory". significativo
cbind.data.frame(model= types_of_model2,
QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_wgt_list[[i]])),
as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_nowgt_list[[i]])),
as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_wgt_list[[i]])),
as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)
),
rbind.data.frame(
dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
print(.)
}
i=3
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>%
{
copiar_nombres2(.)
print(.)
}
cbind.data.frame(model= types_of_model2,
QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_wgt_list[[i]])),
as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_nowgt_list[[i]])),
as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_wgt_list[[i]])),
as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)
),
rbind.data.frame(
dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
print(.)
}
i=4
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>%
{
copiar_nombres2(.)
print(.)
}
cbind.data.frame(model= types_of_model2,
QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_wgt_list[[i]])),
as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_nowgt_list[[i]])),
as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_wgt_list[[i]])),
as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)
),
rbind.data.frame(
dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
print(.)
}
i=5
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>%
{
copiar_nombres2(.)
print(.)
}
cbind.data.frame(model= types_of_model2,
QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_wgt_list[[i]])),
as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_nowgt_list[[i]])),
as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_wgt_list[[i]])),
as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)
),
rbind.data.frame(
dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
))%>%
{
copiar_nombres2(.)
assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
print(.)
}Code
rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::mutate(min_qic_setting=ifelse(QIC==min(QIC),1,0)) %>%
dplyr::mutate(QIC= format(as.numeric(sprintf("%1.1f",QIC)), big.mark=",")) %>%
dplyr::ungroup() %>%
dplyr::mutate(`RR (95% IC)`= paste0(Estimates, " (", `CI.lower`, ", ", CI.upper, ")")) %>%
dplyr::select(setting, model, QIC, `RR (95% IC)`, p.value) %>%
{
#copiar_nombres2()
write.table(., file = paste0(folder_path,"gee_240515.csv"), dec=",", sep="\t")
knitr::kable(., size=10, format="html", caption="Relative risk of treatment non-completion status by reported polysubstance use") %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| setting | model | QIC | RR (95% IC) | p.value |
|---|---|---|---|---|
| basic ambulatory | No time-varying coefficients, no weight | 19,453.6 | 1.03 (1.00, 1.05) | 0.0689 |
| basic ambulatory | No time-varying coefficients, no weight, NB | 19,838.6 | 1.03 (1.00, 1.05) | 0.0698 |
| basic ambulatory | No time-varying coefficients, lag=0 | 19,465.1 | 1.02 (0.99, 1.05) | 0.1648 |
| basic ambulatory | No time-varying coefficients, lag=0, NB | 19,465.1 | 1.02 (0.99, 1.05) | 0.1648 |
| basic ambulatory | No time-varying coefficients, lag=1 | 19,458.8 | 1.02 (1.00, 1.05) | 0.0779 |
| basic ambulatory | No time-varying coefficients, lag=1, NB | 19,458.8 | 1.02 (1.00, 1.05) | 0.0779 |
| basic ambulatory | With time-varying coefficients, no weight | 19,453.6 | 1.03 (1.00, 1.05) | 0.0689 |
| basic ambulatory | With time-varying coefficients, no weight, NB | 19,453.6 | 1.03 (1.00, 1.05) | 0.0689 |
| basic ambulatory | With time-varying coefficients, lag=0 | 19,479.4 | 1.03 (0.99, 1.06) | 0.1328 |
| basic ambulatory | With time-varying coefficients, lag=0, NB | 19,479.4 | 1.03 (0.99, 1.06) | 0.1328 |
| basic ambulatory | With time-varying coefficients, lag=1 | 19,519.4 | 1.03 (0.98, 1.07) | 0.2275 |
| basic ambulatory | With time-varying coefficients, lag=1, NB | 19,519.4 | 1.03 (0.98, 1.07) | 0.2275 |
| basic ambulatory | Stratified by follow-up intervals, no weight | 19,453.6 | 1.03 (1.00, 1.05) | 0.0689 |
| basic ambulatory | Stratified by follow-up intervals, no weight, NB | 19,453.6 | 1.03 (1.00, 1.05) | 0.0689 |
| basic ambulatory | Stratified by follow-up intervals, lag=0 | 19,476.8 | 1.02 (0.99, 1.06) | 0.2682 |
| basic ambulatory | Stratified by follow-up intervals, lag=0, NB | 19,476.8 | 1.02 (0.99, 1.06) | 0.2682 |
| basic ambulatory | Stratified by follow-up intervals, lag=1 | 19,467.2 | 1.04 (1.01, 1.07) | 0.0202 |
| basic ambulatory | Stratified by follow-up intervals, lag=1, NB | 19,467.2 | 1.04 (1.01, 1.07) | 0.0202 |
| basic ambulatory | Stratified by follow-up intervals, 2nd, lag=1 | 19,505.4 | 1.03 (0.99, 1.07) | 0.1745 |
| basic ambulatory | Stratified by follow-up intervals, 2nd, lag=1, NB | 19,505.4 | 1.03 (0.99, 1.07) | 0.1745 |
| GP intensive ambulatory | No time-varying coefficients, no weight | 22,261.9 | 1.04 (1.01, 1.07) | 0.0082 |
| GP intensive ambulatory | No time-varying coefficients, no weight, NB | 22,686.8 | 1.04 (1.01, 1.07) | 0.008 |
| GP intensive ambulatory | No time-varying coefficients, lag=0 | 22,275.9 | 1.04 (1.01, 1.08) | 0.01 |
| GP intensive ambulatory | No time-varying coefficients, lag=0, NB | 22,276.0 | 1.04 (1.01, 1.08) | 0.01 |
| GP intensive ambulatory | No time-varying coefficients, lag=1 | 22,271.9 | 1.04 (1.01, 1.07) | 0.0113 |
| GP intensive ambulatory | No time-varying coefficients, lag=1, NB | 22,271.9 | 1.04 (1.01, 1.07) | 0.0113 |
| GP intensive ambulatory | With time-varying coefficients, no weight | 22,261.9 | 1.04 (1.01, 1.07) | 0.0082 |
| GP intensive ambulatory | With time-varying coefficients, no weight, NB | 22,261.9 | 1.04 (1.01, 1.07) | 0.0082 |
| GP intensive ambulatory | With time-varying coefficients, lag=0 | 22,284.4 | 1.06 (1.02, 1.09) | 0.0032 |
| GP intensive ambulatory | With time-varying coefficients, lag=0, NB | 22,284.4 | 1.06 (1.02, 1.09) | 0.0032 |
| GP intensive ambulatory | With time-varying coefficients, lag=1 | 22,324.6 | 1.05 (1.00, 1.10) | 0.0443 |
| GP intensive ambulatory | With time-varying coefficients, lag=1, NB | 22,324.6 | 1.05 (1.00, 1.10) | 0.0443 |
| GP intensive ambulatory | Stratified by follow-up intervals, no weight | 22,261.9 | 1.04 (1.01, 1.07) | 0.0082 |
| GP intensive ambulatory | Stratified by follow-up intervals, no weight, NB | 22,261.9 | 1.04 (1.01, 1.07) | 0.0082 |
| GP intensive ambulatory | Stratified by follow-up intervals, lag=0 | 22,283.0 | 1.04 (1.01, 1.08) | 0.0174 |
| GP intensive ambulatory | Stratified by follow-up intervals, lag=0, NB | 22,283.0 | 1.04 (1.01, 1.08) | 0.0174 |
| GP intensive ambulatory | Stratified by follow-up intervals, lag=1 | 22,269.4 | 1.02 (0.99, 1.05) | 0.1429 |
| GP intensive ambulatory | Stratified by follow-up intervals, lag=1, NB | 22,269.4 | 1.02 (0.99, 1.05) | 0.1429 |
| GP intensive ambulatory | Stratified by follow-up intervals, 2nd, lag=1 | 22,280.0 | 1.01 (0.98, 1.05) | 0.4077 |
| GP intensive ambulatory | Stratified by follow-up intervals, 2nd, lag=1, NB | 22,280.0 | 1.01 (0.98, 1.05) | 0.4077 |
| GP residential | No time-varying coefficients, no weight | 9,818.1 | 0.97 (0.92, 1.02) | 0.2207 |
| GP residential | No time-varying coefficients, no weight, NB | 9,979.2 | 0.97 (0.92, 1.02) | 0.2203 |
| GP residential | No time-varying coefficients, lag=0 | 9,833.7 | 0.97 (0.92, 1.02) | 0.2773 |
| GP residential | No time-varying coefficients, lag=0, NB | 9,833.7 | 0.97 (0.92, 1.02) | 0.2773 |
| GP residential | No time-varying coefficients, lag=1 | 9,827.8 | 0.95 (0.90, 1.01) | 0.1236 |
| GP residential | No time-varying coefficients, lag=1, NB | 9,827.8 | 0.95 (0.90, 1.01) | 0.1236 |
| GP residential | With time-varying coefficients, no weight | 9,818.1 | 0.97 (0.92, 1.02) | 0.2207 |
| GP residential | With time-varying coefficients, no weight, NB | 9,818.1 | 0.97 (0.92, 1.02) | 0.2207 |
| GP residential | With time-varying coefficients, lag=0 | 9,830.3 | 0.98 (0.93, 1.04) | 0.4692 |
| GP residential | With time-varying coefficients, lag=0, NB | 9,830.3 | 0.98 (0.93, 1.04) | 0.4692 |
| GP residential | With time-varying coefficients, lag=1 | 9,868.4 | 1.00 (0.94, 1.07) | 0.991 |
| GP residential | With time-varying coefficients, lag=1, NB | 9,868.4 | 1.00 (0.94, 1.07) | 0.991 |
| GP residential | Stratified by follow-up intervals, no weight | 9,818.1 | 0.97 (0.92, 1.02) | 0.2207 |
| GP residential | Stratified by follow-up intervals, no weight, NB | 9,818.1 | 0.97 (0.92, 1.02) | 0.2207 |
| GP residential | Stratified by follow-up intervals, lag=0 | 9,839.4 | 0.95 (0.89, 1.02) | 0.137 |
| GP residential | Stratified by follow-up intervals, lag=0, NB | 9,839.4 | 0.95 (0.89, 1.02) | 0.137 |
| GP residential | Stratified by follow-up intervals, lag=1 | 9,837.5 | 0.98 (0.91, 1.04) | 0.4662 |
| GP residential | Stratified by follow-up intervals, lag=1, NB | 9,837.5 | 0.98 (0.91, 1.04) | 0.4662 |
| GP residential | Stratified by follow-up intervals, 2nd, lag=1 | 9,885.3 | 1.01 (0.92, 1.11) | 0.8284 |
| GP residential | Stratified by follow-up intervals, 2nd, lag=1, NB | 9,885.3 | 1.01 (0.92, 1.11) | 0.8284 |
| WO intensive ambulatory | No time-varying coefficients, no weight | 3,424.2 | 0.99 (0.93, 1.05) | 0.715 |
| WO intensive ambulatory | No time-varying coefficients, no weight, NB | 3,487.8 | 0.99 (0.92, 1.06) | 0.7171 |
| WO intensive ambulatory | No time-varying coefficients, lag=0 | 3,434.7 | 0.99 (0.92, 1.07) | 0.8848 |
| WO intensive ambulatory | No time-varying coefficients, lag=0, NB | 3,434.7 | 0.99 (0.92, 1.07) | 0.8848 |
| WO intensive ambulatory | No time-varying coefficients, lag=1 | 3,431.3 | 0.99 (0.92, 1.06) | 0.7415 |
| WO intensive ambulatory | No time-varying coefficients, lag=1, NB | 3,431.3 | 0.99 (0.92, 1.06) | 0.7415 |
| WO intensive ambulatory | With time-varying coefficients, no weight | 3,424.2 | 0.99 (0.93, 1.05) | 0.715 |
| WO intensive ambulatory | With time-varying coefficients, no weight, NB | 3,424.2 | 0.99 (0.93, 1.05) | 0.715 |
| WO intensive ambulatory | With time-varying coefficients, lag=0 | 3,445.9 | 1.01 (0.92, 1.11) | 0.8115 |
| WO intensive ambulatory | With time-varying coefficients, lag=0, NB | 3,445.9 | 1.01 (0.92, 1.11) | 0.8115 |
| WO intensive ambulatory | With time-varying coefficients, lag=1 | 3,468.6 | 1.01 (0.91, 1.12) | 0.8504 |
| WO intensive ambulatory | With time-varying coefficients, lag=1, NB | 3,468.6 | 1.01 (0.91, 1.12) | 0.8504 |
| WO intensive ambulatory | Stratified by follow-up intervals, no weight | 3,424.2 | 0.99 (0.93, 1.05) | 0.715 |
| WO intensive ambulatory | Stratified by follow-up intervals, no weight, NB | 3,424.2 | 0.99 (0.93, 1.05) | 0.715 |
| WO intensive ambulatory | Stratified by follow-up intervals, lag=0 | 3,453.0 | 0.98 (0.90, 1.07) | 0.6456 |
| WO intensive ambulatory | Stratified by follow-up intervals, lag=0, NB | 3,453.0 | 0.98 (0.90, 1.07) | 0.6456 |
| WO intensive ambulatory | Stratified by follow-up intervals, lag=1 | 3,437.3 | 0.98 (0.91, 1.05) | 0.547 |
| WO intensive ambulatory | Stratified by follow-up intervals, lag=1, NB | 3,437.3 | 0.98 (0.91, 1.05) | 0.547 |
| WO intensive ambulatory | Stratified by follow-up intervals, 2nd, lag=1 | 3,466.2 | 0.94 (0.87, 1.03) | 0.1744 |
| WO intensive ambulatory | Stratified by follow-up intervals, 2nd, lag=1, NB | 3,466.2 | 0.94 (0.87, 1.03) | 0.1744 |
| WO residential | No time-varying coefficients, no weight | 4,826.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | No time-varying coefficients, no weight, NB | 4,909.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | No time-varying coefficients, lag=0 | 4,839.1 | 1.15 (1.06, 1.25) | 0.001 |
| WO residential | No time-varying coefficients, lag=0, NB | 4,839.1 | 1.15 (1.06, 1.25) | 0.001 |
| WO residential | No time-varying coefficients, lag=1 | 4,835.6 | 1.13 (1.04, 1.22) | 0.0027 |
| WO residential | No time-varying coefficients, lag=1, NB | 4,835.6 | 1.13 (1.04, 1.22) | 0.0027 |
| WO residential | With time-varying coefficients, no weight | 4,826.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | With time-varying coefficients, no weight, NB | 4,826.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | With time-varying coefficients, lag=0 | 4,839.7 | 1.11 (1.02, 1.21) | 0.0134 |
| WO residential | With time-varying coefficients, lag=0, NB | 4,839.7 | 1.11 (1.02, 1.21) | 0.0134 |
| WO residential | With time-varying coefficients, lag=1 | 4,861.0 | 1.13 (1.03, 1.25) | 0.012 |
| WO residential | With time-varying coefficients, lag=1, NB | 4,861.0 | 1.13 (1.03, 1.25) | 0.012 |
| WO residential | Stratified by follow-up intervals, no weight | 4,826.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | Stratified by follow-up intervals, no weight, NB | 4,826.1 | 1.14 (1.06, 1.23) | 6e-04 |
| WO residential | Stratified by follow-up intervals, lag=0 | 4,851.1 | 1.12 (1.03, 1.22) | 0.0105 |
| WO residential | Stratified by follow-up intervals, lag=0, NB | 4,851.1 | 1.12 (1.03, 1.22) | 0.0105 |
| WO residential | Stratified by follow-up intervals, lag=1 | 4,841.8 | 1.11 (1.02, 1.20) | 0.0117 |
| WO residential | Stratified by follow-up intervals, lag=1, NB | 4,841.8 | 1.11 (1.02, 1.20) | 0.0117 |
| WO residential | Stratified by follow-up intervals, 2nd, lag=1 | 4,874.2 | 1.09 (0.99, 1.20) | 0.0693 |
| WO residential | Stratified by follow-up intervals, 2nd, lag=1, NB | 4,874.2 | 1.09 (0.99, 1.20) | 0.0693 |
Stratifying by presence of alcohol among poly-substance use
Code
data_mine_miss_restr_proc2exp2<-
data_mine_miss_restr_proc2 %>%
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
tidylog::left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","fech_ing_num","otras_sus1_mod","otras_sus2_mod","otras_sus3_mod")], by=c("hash_key"="hash_key","fech_ing_num"="fech_ing_num")) %>%
{if (nrow(.) > nrow(data_mine_miss_restr_proc2)) {message("left join added rows"); .} else .} %>%
dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>%
dplyr::mutate(
susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
) %>%
rowwise() %>%
#Condition
dplyr::mutate(oh_otras_sus = if_else("Alcohol" %in% c(otras_sus1_mod, otras_sus2_mod, otras_sus3_mod),1,0)) %>%
dplyr::mutate(policonsumo2_rec= dplyr::case_when(policonsumo2==1 & oh_otras_sus==1~"2.both",policonsumo2==1 & oh_otras_sus==0~"1.only PSU",T~"0.no PSU")) %>%
ungroup()Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#
nocorr_nowgt_int_list<-list()
nocorr_nowgt_int_probs_list<-list()
nocor_nowgt_int_tab <-list()
for (i in seq_along(plan_names)) {
# Subset the data for the current plan type
current_data <- subset(data_mine_miss_restr_proc2exp2, tipo_de_plan_2_mod == plan_names[i])
# Fit the GEE model for the current subset
model <- geeglm(tr_outcome ~ policonsumo2_rec +
comp_bpsc_y3_severe+
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr, #origen_ingreso_mod fis_comorbidity_icd_10 dg_trs_cons_sus_or
id = id,
data = current_data,
family = poisson,
corstr = "independence")
# Assign the model to the list with a name based on the plan name
model_name <- gsub(" ", "_", plan_names[i])
model_name <- gsub("[^[:alnum:]_]", "", model_name) # Clean up non-alphanumeric characters
nocor_nowgt_int_tab[[paste0("tab_", model_name)]] <- table(current_data$policonsumo2_rec, current_data$tr_outcome) %>% data.frame()
nocorr_nowgt_int_list[[paste("model", model_name, sep = "_")]] <- model
emmeans_response <- emmeans(model, ~policonsumo2_rec, rg.limit=2e5, type = "response")
nocorr_nowgt_int_probs_list[[paste0("emmeans_response_int_", model_name)]] <- emmeans_response
assign(paste0("emmeans_response_int_", model_name), emmeans_response, envir = .GlobalEnv)
}
nocor_nowgt_int_tab[1]$tab_basic_ambulatory
Var1 Var2 Freq
1 0.no PSU 0 629
2 1.only PSU 0 1210
3 2.both 0 286
4 0.no PSU 1 1827
5 1.only PSU 1 4450
6 2.both 1 1591Code
nocor_nowgt_int_tab[2]$tab_GP_intensive_ambulatory
Var1 Var2 Freq
1 0.no PSU 0 624
2 1.only PSU 0 1615
3 2.both 0 401
4 0.no PSU 1 1573
5 1.only PSU 1 5226
6 2.both 1 2058Code
nocor_nowgt_int_tab[3]$tab_GP_residential
Var1 Var2 Freq
1 0.no PSU 0 261
2 1.only PSU 0 905
3 2.both 0 400
4 0.no PSU 1 572
5 1.only PSU 1 2360
6 2.both 1 706Code
nocor_nowgt_int_tab[4]$tab_WO_intensive_ambulatory
Var1 Var2 Freq
1 0.no PSU 0 110
2 1.only PSU 0 253
3 2.both 0 61
4 0.no PSU 1 292
5 1.only PSU 1 726
6 2.both 1 318Code
nocor_nowgt_int_tab[5]$tab_WO_residential
Var1 Var2 Freq
1 0.no PSU 0 172
2 1.only PSU 0 413
3 2.both 0 130
4 0.no PSU 1 289
5 1.only PSU 1 1147
6 2.both 1 383Code
# > nocor_nowgt_int_tab[1]
# $tab_basic_ambulatory
# Var1 Var2 Freq
# 1 0.no PSU 0 629
# 2 1.only PSU 0 1210
# 3 2.both 0 286
# 4 0.no PSU 1 1827
# 5 1.only PSU 1 4450
# 6 2.both 1 1591
#
# > nocor_nowgt_int_tab[2]
# $tab_GP_intensive_ambulatory
# Var1 Var2 Freq
# 1 0.no PSU 0 624
# 2 1.only PSU 0 1615
# 3 2.both 0 401
# 4 0.no PSU 1 1573
# 5 1.only PSU 1 5226
# 6 2.both 1 2058
#
# > nocor_nowgt_int_tab[3]
# $tab_GP_residential
# Var1 Var2 Freq
# 1 0.no PSU 0 261
# 2 1.only PSU 0 905
# 3 2.both 0 400
# 4 0.no PSU 1 572
# 5 1.only PSU 1 2360
# 6 2.both 1 706
#
# > nocor_nowgt_int_tab[4]
# $tab_WO_intensive_ambulatory
# Var1 Var2 Freq
# 1 0.no PSU 0 110
# 2 1.only PSU 0 253
# 3 2.both 0 61
# 4 0.no PSU 1 292
# 5 1.only PSU 1 726
# 6 2.both 1 318
#
# > nocor_nowgt_int_tab[5]
# $tab_WO_residential
# Var1 Var2 Freq
# 1 0.no PSU 0 172
# 2 1.only PSU 0 413
# 3 2.both 0 130
# 4 0.no PSU 1 289
# 5 1.only PSU 1 1147
# 6 2.both 1 383
#
invisible("no hay tan pocos casos en un estrato, salvo en WO intensive ambulatory")
invisible("En el basic ambulatory, GP intensive-amb y GP residential, el sig es el ambos (consumo alcohol)")
#broom::tidy(nocorr_nowgt_int_list[[1]], exponentiate=T, conf.int=T)[1:3,]
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 0.00000870 4.83 5.81 0.0159 6.68e-10 0.113
# 2 policonsumo2_rec1.only PSU 1.01 0.0143 0.488 0.485 9.82e- 1 1.04
# 3 policonsumo2_rec2.both 1.08 0.0161 24.8 0.000000622 1.05e+ 0 1.12
#broom::tidy(nocorr_nowgt_int_list[[2]], exponentiate=T, conf.int=T)[1:3,]
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 0.00000655 4.81 6.15 0.0132 5.23e-10 0.0822
# 2 policonsumo2_rec1.only PSU 1.02 0.0153 1.98 0.160 9.92e- 1 1.05
# 3 policonsumo2_rec2.both 1.10 0.0162 35.1 0.00000000312 1.07e+ 0 1.14
invisible("Aunque GP residential es protector el policonsumo en both")
#broom::tidy(nocorr_nowgt_int_list[[3]], exponentiate=T, conf.int=T)[1:3,]
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 1.14 9.06 0.000199 0.989 0.0000000221 58484383.
# 2 policonsumo2_rec1.only PSU 0.989 0.0275 0.175 0.676 0.937 1.04
# 3 policonsumo2_rec2.both 0.885 0.0332 13.5 0.000243 0.830 0.945
invisible("En WO intensive ambulatory, no tira ninguno")
#broom::tidy(nocorr_nowgt_int_list[[4]], exponentiate=T, conf.int=T)[1:3,]
# # A tibble: 3 x 7
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 0.000318 12.0 0.452 0.501 2.02e-14 4991312.
# 2 policonsumo2_rec1.only PSU 0.956 0.0357 1.62 0.203 8.91e- 1 1.02
# 3 policonsumo2_rec2.both 1.07 0.0376 3.13 0.0767 9.93e- 1 1.15
#
invisible("En WO residential, ambos mueven la aguja a mayor riesgo")
#
#broom::tidy(nocorr_nowgt_int_list[[5]], exponentiate=T, conf.int=T)[1:3,]
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 249351. 11.6 1.15 0.284 0.0000331 1.88e15
# 2 policonsumo2_rec1.only PSU 1.14 0.0388 11.3 0.000775 1.06 1.23e 0
# 3 policonsumo2_rec2.both 1.14 0.0439 8.83 0.00297 1.05 1.24e 0
#
#
bind_rows(
lapply(nocorr_nowgt_int_list[1:5], function(model) {
broom::tidy(model, exponentiate = TRUE, conf.int = TRUE) %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value) %>%
slice(2:3)
}),
.id = "model"
) %>%
dplyr::mutate(term= gsub("policonsumo2_","", term)) %>%
dplyr::mutate(across(c("estimate","conf.low", "conf.high"),~sprintf("%1.2f",.))) %>%
dplyr::mutate(RR= paste0(estimate, " (", conf.low,"-", conf.high,")")) %>%
dplyr::select(-estimate, -conf.low, -conf.high) %>%
dplyr::select(model, term, RR, p.value) %>%
dplyr::mutate(across(c("p.value"),~sprintf("%1.4f",.))) %>%
{
knitr::kable(., "markdown", caption= "Relative risk of treatment non-completion status by reported polysubstance use, whether ")
}| model | term | RR | p.value |
|---|---|---|---|
| model_basic_ambulatory | rec1.only PSU | 1.01 (0.98-1.04) | 0.4849 |
| model_basic_ambulatory | rec2.both | 1.08 (1.05-1.12) | 0.0000 |
| model_GP_intensive_ambulatory | rec1.only PSU | 1.02 (0.99-1.05) | 0.1599 |
| model_GP_intensive_ambulatory | rec2.both | 1.10 (1.07-1.14) | 0.0000 |
| model_GP_residential | rec1.only PSU | 0.99 (0.94-1.04) | 0.6761 |
| model_GP_residential | rec2.both | 0.89 (0.83-0.94) | 0.0002 |
| model_WO_intensive_ambulatory | rec1.only PSU | 0.96 (0.89-1.02) | 0.2034 |
| model_WO_intensive_ambulatory | rec2.both | 1.07 (0.99-1.15) | 0.0767 |
| model_WO_residential | rec1.only PSU | 1.14 (1.06-1.23) | 0.0008 |
| model_WO_residential | rec2.both | 1.14 (1.05-1.24) | 0.0030 |
Code
# bind_rows(
# lapply(1:5, function(model) {
# nocor_nowgt_int_tab[[model]]%>% dplyr::mutate(Var0=ifelse(grepl("no PSU",Var1),1,0)) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n= sum(Freq))
# }),
# .id = "model"
# ) %>%
# group_by(model) %>%
# summarise(n=sum(n))
data_barplot <- bind_rows(
lapply(1:5, function(model) {
nocor_nowgt_int_tab[[model]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq))
}),
.id = "model"
) %>%
group_by(model, Var0) %>%
mutate(percentage = n / sum(n)) %>%
ungroup()Code
# Filter data for treatment non-completion only
data_non_completion <- data_barplot %>% filter(Var2 == "1")
# Create labels for the total counts and percentages for each group
labels <- data_barplot %>%
group_by(model, Var0) %>%
summarise(total = sum(n)) %>%
mutate(label = scales::percent(total / sum(total), accuracy = 0.1))Code
# Calculate overall percentages for PSU and no PSU
overall_percentages <- data_barplot %>%
group_by(model, Var0) %>%
summarise(total = sum(n)) %>%
mutate(overall_percentage = total / sum(total))Code
# Create a separate data frame for the labels
labels <- overall_percentages %>%
mutate(label = scales::percent(overall_percentage, accuracy = 1))
data_barplot$Var0 <- factor(data_barplot$Var0, levels = c("Single\nsubstance use", "Polysubstance\nuse"))
labels$Var0 <- factor(labels$Var0, levels = c("Single\nsubstance use", "Polysubstance\nuse"))
showtext_auto()
showtext_opts(dpi = 500)
invisible("Si no funciona el times new roman")
#windowsFonts(TimesNewRoman = windowsFont("Times New Roman"))
# names(grDevices::windowsFonts())
# Generate the plot
ggplot(data_barplot, aes(x = factor(Var0), y = n, fill = factor(Var2))) +
geom_bar(stat = "identity", position = position_stack(reverse = TRUE)) +
#geom_bar(stat = "identity", position = "stack") +
facet_wrap(~ model, ncol = 1, labeller = labeller(model = c(
"1" = "Basic ambulatory (n= 9,993)",
"2" = "General-population, intensive ambulatory (n= 11,497)",
"3" = "General-population, residential (n= 5,204)",
"4" = "Women-specific, intensive ambulatory (n= 1,760)",
"5" = "Women-specific, residential (n= 2,534)"
)), strip.position = "top") +
scale_fill_manual(values = c("grey80", "grey35"), label=c("Completion", "Non-completion")) +
labs(x = "Polysubstance use (PSU)", y = "Count", fill = "Treatment\noutcome") +
geom_label(data = data_non_completion, aes(label = scales::percent(percentage, accuracy = 1)), position = position_stack(vjust = 0.8, reverse = TRUE), size = 3.9, fontface = "bold", color = "white", fill = "black", label.size = 0.2, family = "serif") +
geom_text(data = labels, aes(x = Var0, y = 9e3, label = label, fill="1"),
size = 3.9, color = "black", hjust = 0.5, vjust = -0.5, fontface = "bold.italic", family = "serif") +
ylim(0,1.1e4)+
#theme_minimal(base_size = 14, base_family = "Times New Roman") +
theme(
strip.text = element_text(size = 12, face = "bold", family = "serif"),
#axis.text.x = element_text(angle = 0, hjust = 1),
axis.title.x = element_blank(),
axis.text = element_text(size = 12, family = "serif"),
axis.text.x = element_text(family = "serif"),
text = element_text(family = "serif"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
strip.text.x = element_text(angle = 0, hjust = 0.5, size = 12, face = "bold", family = "serif"),
strip.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
legend.text = element_text(size= 12, family = "serif"),
legend.title = element_text(size= 13, family = "serif"),
)Warning in geom_text(data = labels, aes(x = Var0, y = 9000, label = label, : Ignoring unknown aesthetics: fill
Code
ggsave("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/_figs/barplot3.png", height=14*.68, width=10*.68, dpi=500)
ggsave("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/_figs/barplot3.jpg", height=14*.68, width=10*.68, dpi=500)
Code
#data_barplot
#
mantelhaen.test(
xtabs(n ~Var2+ Var0 +model, data= data_barplot)
)Warning in (x[1L, 1L, ] + x[2L, 2L, ]) * x[1L, 1L, ] * x[2L, 2L, ]: NAs producidos por enteros excedidos
Warning in (x[1L, 1L, ] + x[2L, 2L, ]) * x[1L, 2L, ] * x[2L, 1L, ]: NAs producidos por enteros excedidos
Warning in (x[1L, 2L, ] + x[2L, 1L, ]) * x[1L, 1L, ] * x[2L, 2L, ]: NAs producidos por enteros excedidos
Warning in (x[1L, 2L, ] + x[2L, 1L, ]) * x[1L, 2L, ] * x[2L, 1L, ]: NAs producidos por enteros excedidos
Mantel-Haenszel chi-squared test with continuity correction
data: xtabs(n ~ Var2 + Var0 + model, data = data_barplot)
Mantel-Haenszel X-squared = 93.675, df = 1, p-value < 2.2e-16
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
NA NA
sample estimates:
common odds ratio
1.361499 Code
#Mantel-Haenszel X-squared = 93.675, df = 1, p-value < 2.2e-16
invisible("Si bien se ve una aso")
chisq.test(xtabs(n ~Var2+ Var0, data= data_barplot))
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = data_barplot)
X-squared = 76.039, df = 1, p-value < 2.2e-16Code
#X-squared = 76.039, df = 1, p-value < 2.2e-16
vcd::woolf_test(
xtabs(n ~Var2+ Var0 +model, data= data_barplot)
)
Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
data: xtabs(n ~ Var2 + Var0 + model, data = data_barplot)
X-squared = 13.735, df = 4, p-value = 0.008192Code
#Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
#X-squared = 13.735, df = 4, p-value = 0.008192
invisible("Los odds de no completar varían por baseline tr. setting (p<0,001)")
invisible("No se puede usar Maentel-haenzel")
invisible("Suggests conditional dependence and varying strength or direction of association across strata (interaction), cautioning against a simple summary of the association.[common OR may not be reliable]")
# nocor_nowgt_int_tab[[1]] %>%
# dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
# dplyr::group_by(Var0, Var2) %>%
# dplyr::summarise(n = sum(Freq))
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[1]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq)))
)
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[1]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 36.389, df = 1, p-value = 1.616e-09Code
#X-squared = 36.389, df = 1, p-value = 1.616e-09
# Phi-Coefficient : 0.061
# Contingency Coeff.: 0.061
# Cramer's V : 0.061
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[2]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq)))
)
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[2]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 45.055, df = 1, p-value = 1.916e-11Code
#X-squared = 45.055, df = 1, p-value = 1.916e-11
# Phi-Coefficient : 0.063
# Contingency Coeff.: 0.063
# Cramer's V : 0.063
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[3]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq)))
)
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[3]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 0.65673, df = 1, p-value = 0.4177Code
#X-squared = 0.65673, df = 1, p-value = 0.4177
# hi-Coefficient : 0.012
# Contingency Coeff.: 0.012
# Cramer's V : 0.012
#
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[4]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq)))
)
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[4]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 2.8231, df = 1, p-value = 0.09291Code
#X-squared = 2.8231, df = 1, p-value = 0.09291
# Phi-Coefficient : 0.042
# Contingency Coeff.: 0.042
# Cramer's V : 0.042
#
invisible("setting 5= No completar vs. completar, según PSU vs. no-PSU")
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[5]] %>%
dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%
dplyr::group_by(Var0, Var2) %>%
dplyr::summarise(n = sum(Freq)))
)
Pearson's Chi-squared test with Yates' continuity correction
data: xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[5]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 22.463, df = 1, p-value = 2.142e-06Code
#X-squared = 22.463, df = 1, p-value = 2.142e-06
# Phi-Coefficient : 0.095
# Contingency Coeff.: 0.095
# Cramer's V : 0.095
invisible("compare setting vs. PSUs")
chisq.test(xtabs(n ~Var0+ model, data= data_barplot))
Pearson's Chi-squared test
data: xtabs(n ~ Var0 + model, data = data_barplot)
X-squared = 194.31, df = 4, p-value < 2.2e-16Code
#X-squared = 194.31, df = 4, p-value < 2.2e-16
# Cox, M. K., & Key, C. H. (1993). Post hoc pair-wise comparisons for the chi-square test of homogeneity of proportions. Educational and Psychological Measurement, 53(4), 951–962. https://doi.org/10.1177/0013164493053004008
invisible("Settings vs. noncompletion")
chisq.test(xtabs(n ~Var2+ model, data= data_barplot))
Pearson's Chi-squared test
data: xtabs(n ~ Var2 + model, data = data_barplot)
X-squared = 177.64, df = 4, p-value < 2.2e-16Code
#X-squared = 177.64, df = 4, p-value < 2.2e-16
cat("For each baseline setting, calculate the std. residuals of chi-square, to compare non-completion shares")For each baseline setting, calculate the std. residuals of chi-square, to compare non-completion sharesCode
data.frame(chisq.posthoc.test(xtabs(n ~Var2+ model, data= data_barplot))) |>
dplyr::mutate(X4=as.character(X4)) |>
pivot_longer(cols = c(`X1`,`X2`,`X3`,`X4`,`X5`), names_to = "Group", values_to = "Values") %>%
pivot_wider(names_from = Value, values_from = Values) %>%
dplyr::filter(Dimension=="1") %>%
dplyr::mutate(across(1, ~ gsub("\n"," ",as.character(.))))|>
dplyr::mutate(Residuals= sprintf("%1.2f", as.numeric(Residuals))) %>%
dplyr::mutate( `p values`= sprintf("%1.3f", as.numeric(gsub("\\*","",`p values`)))) %>%
dplyr::select(Dimension, Residuals, `p values`) %>%
dplyr::mutate(`p values`= ifelse(`p values`=="0.000","<0.001",`p values`))# A tibble: 5 x 3
Dimension Residuals `p values`
<chr> <chr> <chr>
1 1 8.07 <0.001
2 1 3.61 0.003
3 1 -11.07 <0.001
4 1 0.02 1.000
5 1 -5.05 <0.001 Code
data.frame(chisq.posthoc.test(xtabs(n ~Var0+ model, data= data_barplot))) |>
dplyr::mutate(X4=as.character(X4)) |>
pivot_longer(cols = c(`X1`,`X2`,`X3`,`X4`,`X5`), names_to = "Group", values_to = "Values") %>%
pivot_wider(names_from = Value, values_from = Values) %>%
dplyr::filter(Dimension=="Polysubstance\nuse") %>%
dplyr::mutate(across(1, ~ gsub("\n"," ",as.character(.))))|>
dplyr::mutate(Residuals= sprintf("%1.2f", as.numeric(Residuals))) %>%
dplyr::mutate( `p values`= sprintf("%1.3f", as.numeric(gsub("\\*","",`p values`)))) %>%
dplyr::select(Dimension, Residuals, `p values`) %>%
dplyr::mutate(`p values`= ifelse(`p values`=="0.000","<0.001",`p values`))# A tibble: 5 x 3
Dimension Residuals `p values`
<chr> <chr> <chr>
1 Polysubstance use -12.30 <0.001
2 Polysubstance use 4.62 <0.001
3 Polysubstance use 8.78 <0.001
4 Polysubstance use -2.52 0.118
5 Polysubstance use 2.99 0.028 Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#
nocorr_nowgt_int2_list<-
bind_rows(
lapply(nocorr_nowgt_int_list[1:5], function(model) {
broom::tidy(model, exponentiate = TRUE, conf.int = TRUE) %>%
dplyr::select(term, estimate, std.error, conf.low, conf.high, p.value) %>%
dplyr::mutate(log_rr= log(estimate), log_lcl= log(conf.low), log_ucl= log(conf.high), variance= ((log_ucl-log_rr)/qnorm(0.975))^2)%>%
dplyr::slice(2:3)
}),
.id = "model"
)
log_rr_int<- nocorr_nowgt_int2_list$log_rr
variances_int<- nocorr_nowgt_int2_list$variance
results_int <- expand.grid(i=1:5, j=1:5) %>%
filter(i < j) %>%
mutate(
drrr = log_rr_int[i] - log_rr_int[j],
var_drrr = variances_int[i] + variances_int[j],
se_drrr = sqrt(var_drrr),
z_value = drrr / se_drrr,
p_value = 2 * (1 - pnorm(abs(z_value))),
ratio_estimates = exp(drrr),
lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
) %>%
dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2],
i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]),
j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2],
j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))Marginal probabilities
Loop over each unique plan name to fit a model. We used an interaction of PSU and treatment setting to capture the heterogeneity of PSU by treatment setting.
Code
df_marginal_plots<-
cbind.data.frame(wgt= c("wgt", "iiw_nocorr_st", "iiw_nocorr_alt_st", "wgt", "iiw_after_ph_st", "iiw_after_ph_st"),
df= c("data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp_iiw_after_ph","data_mine_miss_restr_proc2exp_iiw_after_ph","data_mine_miss_restr_proc2exp_iiw_after_ph_alt"),
model= c("no corr no wgt","no corr main","no corr alt", "after ph no wgt", "after ph main", "after ph alt"))
# Loop over each unique plan name to fit a model
for (i in seq_along(df_marginal_plots$wgt)) {
# Subset the data for the current plan type
#current_data <- subset(data_mine_miss_restr_proc2, tipo_de_plan_2_mod == plan_names[i])
geeglm( tr_outcome ~ policonsumo2 * tipo_de_plan_2_mod +
edad_al_ing_1 +
ano_nac_corr +
susinidumrec_otr +
susinidumrec_coc +
susinidumrec_pbc +
susinidumrec_mar +
psycom_dum_study +
psycom_dum_with +
freq_cons_dum_5day +
cond_oc_dum_2inact +
cond_oc_dum_3unemp +
susprindumrec_coc +
susprindumrec_pbc +
susprindumrec_mar +
susprindumrec_otr,
id= id,
data= get(df_marginal_plots$df[i]) %>% dplyr::mutate(wgt=1),
weight= get(df_marginal_plots$df[i]) %>% dplyr::mutate(wgt=1) %>% pull(df_marginal_plots$wgt[i]),
family= poisson,
corstr= "independence") %>%
assign(paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i])), ., envir = .GlobalEnv)
print(message(paste0("Models: ",paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i])))))
emmeans_response <- emmeans(get(paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i]))), ~policonsumo2+tipo_de_plan_2_mod, rg.limit=2e5, type = "response")%>% assign(paste0("emmeans_response_", gsub(" ","_",df_marginal_plots$model[i])), ., envir = .GlobalEnv)
print(message(paste0("Emmeans: ",paste0("emmeans_response_", gsub(" ","_",df_marginal_plots$model[i])))))
}NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULLCode
plot(emmeans_response_no_corr_no_wgt)+ ggtitle("No corr, no weights")#+ geom_point(aes(x = percent, y = conc), data = pigs, pch = 2, color = "blue")
Code
plot(emmeans_response_no_corr_main)+ ggtitle("No corr, main weights")
Code
plot(emmeans_response_no_corr_alt)+ ggtitle("No corr, alt weights")
Code
plot(emmeans_response_after_ph_no_wgt)+ ggtitle("After PH, no weights")
Code
plot(emmeans_response_after_ph_main)+ ggtitle("After PH, main weights")
Code
plot(emmeans_response_after_ph_alt)+ ggtitle("After PH, alt weights")
Heterogeneity
Code
library(meta)
try(library(dmetar))Error in library(dmetar) : there is no package called 'dmetar'Code
library(metafor)
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##
# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`Estimates`) %>% as.numeric()
lcl <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`CI.lower`) %>% as.numeric()
ucl <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`CI.upper`) %>% as.numeric()
# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr <- log(rrr)
log_lcl <- log(lcl)
log_ucl <- log(ucl)
# Calcula las varianzas
variances <- ((log_ucl - log_rrr) / qnorm(0.975))^2
# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis <- rma(yi=log_rrr, sei=sqrt(variances), method="REML")
# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q <- meta_analysis$QE
p_value_q <- meta_analysis$QEp
cat(sprintf("No time-varying coefficients, no weight- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q, p_value_q))No time-varying coefficients, no weight- Cochran's Q = 14.49, p-value = 0.0059Code
#concuerda con el calculo de stata
#
#
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###
###
# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`Estimates`) %>% as.numeric()
lcl2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`CI.lower`) %>% as.numeric()
ucl2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`CI.upper`) %>% as.numeric()
# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr2 <- log(rrr2)
log_lcl2 <- log(lcl2)
log_ucl2 <- log(ucl2)
# Calcula las varianzas
variances2 <- ((log_ucl2 - log_rrr2) / qnorm(0.975))^2
# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis2 <- rma(yi=log_rrr2, sei=sqrt(variances2), method="REML")
# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q2 <- meta_analysis2$QE
p_value_q2 <- meta_analysis2$QEp
cat(sprintf("No time-varying coefficients, lag=0, main- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q2, p_value_q2))No time-varying coefficients, lag=0, main- Cochran's Q = 13.32, p-value = 0.0098Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`Estimates`) %>% as.numeric()
lcl3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`CI.lower`) %>% as.numeric()
ucl3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>%
dplyr::group_by(setting) %>%
dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`CI.upper`) %>% as.numeric()
# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr3 <- log(rrr3)
log_lcl3 <- log(lcl3)
log_ucl3 <- log(ucl3)
# Calcula las varianzas
variances3 <- ((log_ucl3 - log_rrr3) / qnorm(0.975))^2
# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis3 <- rma(yi=log_rrr3, sei=sqrt(variances3), method="REML")
# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q3 <- meta_analysis3$QE
p_value_q3 <- meta_analysis3$QEp
cat(sprintf("No time-varying coefficients, lag=1- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q3, p_value_q3))No time-varying coefficients, lag=1- Cochran's Q = 14.24, p-value = 0.0066Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###
# Comparaciones de Altman para todas las combinaciones de RRR
results <- expand.grid(i=1:5, j=1:5) %>%
filter(i < j) %>%
mutate(
drrr = log_rrr[i] - log_rrr[j],
var_drrr = variances[i] + variances[j],
se_drrr = sqrt(var_drrr),
z_value = drrr / se_drrr,
p_value = 2 * (1 - pnorm(abs(z_value))),
ratio_estimates = exp(drrr),
lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
) %>%
dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))
results2 <- expand.grid(i=1:5, j=1:5) %>%
filter(i < j) %>%
mutate(
drrr = log_rrr2[i] - log_rrr2[j],
var_drrr = variances2[i] + variances2[j],
se_drrr = sqrt(var_drrr),
z_value = drrr / se_drrr,
p_value = 2 * (1 - pnorm(abs(z_value))),
ratio_estimates = exp(drrr),
lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
) %>%
dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))
results3 <- expand.grid(i=1:5, j=1:5) %>%
filter(i < j) %>%
mutate(
drrr = log_rrr3[i] - log_rrr3[j],
var_drrr = variances3[i] + variances3[j],
se_drrr = sqrt(var_drrr),
z_value = drrr / se_drrr,
p_value = 2 * (1 - pnorm(abs(z_value))),
ratio_estimates = exp(drrr),
lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
) %>%
dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))
# Imprime los resultados
#el 5 tiene el mayor efecto en omparación con el resto
#
bind_rows(cbind.data.frame(model= "No time-varying coefficients, no weight",results),
cbind.data.frame(model= "No time-varying coefficients, lag=0",results2),
cbind.data.frame(model= "No time-varying coefficients, lag=1",results3)
) %>%
dplyr::mutate(across(c("drrr","var_drrr", "se_drrr", "z_value", "ratio_estimates", "lower_bound", "upper_bound"),~sprintf("%1.2f",.))) %>%
dplyr::mutate(across(c("p_value"),~sprintf("%1.4f",.))) %>%
dplyr::mutate(`rrr (95% IC)`= paste0(ratio_estimates, " (", lower_bound, ", ", upper_bound,")")) %>%
dplyr::select(-any_of(c("drrr", "var_drrr", "se_drrr", "z_value", "ratio_estimates", "lower_bound", "upper_bound"))) %>%
{
#copiar_nombres2()
write.table(., file = paste0(folder_path,"rr_contrasts_240625.csv"), dec=",", sep="\t")
knitr::kable(., size=10, format="html", caption="GEE Models, Heterogeneity & contrasts by setting") %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px")
}| model | i | j | p_value | rrr (95% IC) |
|---|---|---|---|---|
| No time-varying coefficients, no weight | basic ambulatory | GP intensive ambulatory | 0.5812 | 0.99 (0.96, 1.02) |
| No time-varying coefficients, no weight | basic ambulatory | GP residential | 0.0288 | 1.06 (1.01, 1.12) |
| No time-varying coefficients, no weight | GP intensive ambulatory | GP residential | 0.0180 | 1.07 (1.01, 1.14) |
| No time-varying coefficients, no weight | basic ambulatory | WO intensive ambulatory | 0.2098 | 1.04 (0.98, 1.11) |
| No time-varying coefficients, no weight | GP intensive ambulatory | WO intensive ambulatory | 0.1395 | 1.05 (0.98, 1.12) |
| No time-varying coefficients, no weight | GP residential | WO intensive ambulatory | 0.6052 | 0.98 (0.91, 1.06) |
| No time-varying coefficients, no weight | basic ambulatory | WO residential | 0.0112 | 0.90 (0.84, 0.98) |
| No time-varying coefficients, no weight | GP intensive ambulatory | WO residential | 0.0266 | 0.91 (0.84, 0.99) |
| No time-varying coefficients, no weight | GP residential | WO residential | 0.0005 | 0.85 (0.78, 0.93) |
| No time-varying coefficients, no weight | WO intensive ambulatory | WO residential | 0.0040 | 0.87 (0.79, 0.96) |
| No time-varying coefficients, lag=0 | basic ambulatory | GP intensive ambulatory | 0.4239 | 0.98 (0.94, 1.03) |
| No time-varying coefficients, lag=0 | basic ambulatory | GP residential | 0.0895 | 1.05 (0.99, 1.11) |
| No time-varying coefficients, lag=0 | GP intensive ambulatory | GP residential | 0.0298 | 1.07 (1.01, 1.14) |
| No time-varying coefficients, lag=0 | basic ambulatory | WO intensive ambulatory | 0.4805 | 1.03 (0.95, 1.12) |
| No time-varying coefficients, lag=0 | GP intensive ambulatory | WO intensive ambulatory | 0.2636 | 1.05 (0.96, 1.15) |
| No time-varying coefficients, lag=0 | GP residential | WO intensive ambulatory | 0.6656 | 0.98 (0.89, 1.07) |
| No time-varying coefficients, lag=0 | basic ambulatory | WO residential | 0.0077 | 0.89 (0.81, 0.97) |
| No time-varying coefficients, lag=0 | GP intensive ambulatory | WO residential | 0.0313 | 0.90 (0.83, 0.99) |
| No time-varying coefficients, lag=0 | GP residential | WO residential | 0.0006 | 0.84 (0.77, 0.93) |
| No time-varying coefficients, lag=0 | WO intensive ambulatory | WO residential | 0.0100 | 0.86 (0.77, 0.96) |
| No time-varying coefficients, lag=1 | basic ambulatory | GP intensive ambulatory | 0.3486 | 0.98 (0.94, 1.02) |
| No time-varying coefficients, lag=1 | basic ambulatory | GP residential | 0.0397 | 1.07 (1.00, 1.15) |
| No time-varying coefficients, lag=1 | GP intensive ambulatory | GP residential | 0.0086 | 1.09 (1.02, 1.17) |
| No time-varying coefficients, lag=1 | basic ambulatory | WO intensive ambulatory | 0.4305 | 1.03 (0.96, 1.11) |
| No time-varying coefficients, lag=1 | GP intensive ambulatory | WO intensive ambulatory | 0.1919 | 1.05 (0.98, 1.13) |
| No time-varying coefficients, lag=1 | GP residential | WO intensive ambulatory | 0.3783 | 0.96 (0.88, 1.05) |
| No time-varying coefficients, lag=1 | basic ambulatory | WO residential | 0.0143 | 0.90 (0.83, 0.98) |
| No time-varying coefficients, lag=1 | GP intensive ambulatory | WO residential | 0.0466 | 0.92 (0.85, 1.00) |
| No time-varying coefficients, lag=1 | GP residential | WO residential | 0.0005 | 0.84 (0.76, 0.93) |
| No time-varying coefficients, lag=1 | WO intensive ambulatory | WO residential | 0.0116 | 0.88 (0.79, 0.97) |
Session info
Code
message(paste0("R library: ", Sys.getenv("R_LIBS_USER")))Code
message(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))Code
message(paste0("Editor context: ", path))Code
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
) %>%
knitr::kable(caption = "R packages", format = "html",
col.names = c("Row number", "Package", "Version"),
row.names = FALSE,
align = c("c", "l", "r")) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px") | Row number | Package | Version |
|---|---|---|
| assertthat | 0.2.1 | CRAN (R 4.1.2) |
| backports | 1.4.1 | CRAN (R 4.1.2) |
| base64enc | 0.1-3 | CRAN (R 4.1.1) |
| BiocManager | 1.30.18 | CRAN (R 4.1.3) |
| biostat3 | 0.1.8 | CRAN (R 4.1.3) |
| boot | 1.3-28.1 | CRAN (R 4.1.3) |
| broom | 1.0.1 | CRAN (R 4.1.3) |
| cachem | 1.0.6 | CRAN (R 4.1.2) |
| callr | 3.7.2 | CRAN (R 4.1.3) |
| cellranger | 1.1.0 | CRAN (R 4.1.2) |
| checkmate | 2.1.0 | CRAN (R 4.1.3) |
| chisq.posthoc.test | 0.1.3 | Github (ebbertd/chisq.posthoc.test@186d2ca6bbdba9fc19601aff4696ae1b85e7e0b0) |
| class | 7.3-20 | CRAN (R 4.1.3) |
| cli | 3.4.1 | CRAN (R 4.1.3) |
| clipr | 0.8.0 | CRAN (R 4.1.2) |
| clisymbols | 1.2.0 | CRAN (R 4.1.3) |
| cluster | 2.1.4 | CRAN (R 4.1.3) |
| coda | 0.19-4 | CRAN (R 4.1.2) |
| codetools | 0.2-19 | CRAN (R 4.1.3) |
| colorspace | 2.0-3 | CRAN (R 4.1.2) |
| CompQuadForm | 1.4.3 | CRAN (R 4.1.1) |
| cowplot | 1.1.1 | CRAN (R 4.1.2) |
| crayon | 1.5.2 | CRAN (R 4.1.3) |
| data.table | 1.14.2 | CRAN (R 4.1.2) |
| data.tree | 1.0.0 | CRAN (R 4.1.3) |
| DBI | 1.1.3 | CRAN (R 4.1.3) |
| dbplyr | 2.2.1 | CRAN (R 4.1.3) |
| deldir | 1.0-6 | CRAN (R 4.1.1) |
| DescTools | 0.99.48 | CRAN (R 4.1.3) |
| devtools | 2.4.4 | CRAN (R 4.1.3) |
| DiagrammeR | 1.0.9 | CRAN (R 4.1.3) |
| digest | 0.6.29 | CRAN (R 4.1.3) |
| dplyr | 1.0.10 | CRAN (R 4.1.3) |
| e1071 | 1.7-11 | CRAN (R 4.1.3) |
| easyalluvial | 0.3.2 | CRAN (R 4.1.2) |
| ellipsis | 0.3.2 | CRAN (R 4.1.2) |
| emmeans | 1.8.1-1 | CRAN (R 4.1.3) |
| estimability | 1.4.1 | CRAN (R 4.1.3) |
| evaluate | 0.16 | CRAN (R 4.1.3) |
| Exact | 3.2 | CRAN (R 4.1.3) |
| expm | 1.0-0 | CRAN (R 4.1.2) |
| extrafont | 0.18 | CRAN (R 4.1.3) |
| extrafontdb | 1.0 | CRAN (R 4.1.1) |
| fansi | 1.0.3 | CRAN (R 4.1.3) |
| farver | 2.1.1 | CRAN (R 4.1.3) |
| fastmap | 1.1.0 | CRAN (R 4.1.3) |
| forcats | 0.5.2 | CRAN (R 4.1.3) |
| foreign | 0.8-83 | CRAN (R 4.1.3) |
| Formula | 1.2-4 | CRAN (R 4.1.1) |
| fs | 1.5.2 | CRAN (R 4.1.3) |
| future | 1.28.0 | CRAN (R 4.1.3) |
| future.apply | 1.9.1 | CRAN (R 4.1.3) |
| gargle | 1.2.1 | CRAN (R 4.1.3) |
| geeasy | 0.1.2 | CRAN (R 4.1.2) |
| geeM | 0.10.1 | CRAN (R 4.1.3) |
| geepack | 1.3.9 | CRAN (R 4.1.3) |
| generics | 0.1.3 | CRAN (R 4.1.3) |
| ggalluvial | 0.12.5 | CRAN (R 4.1.3) |
| ggforce | 0.4.1 | CRAN (R 4.1.3) |
| ggformula | 0.10.2 | CRAN (R 4.1.3) |
| ggplot2 | 3.4.1 | CRAN (R 4.1.3) |
| ggridges | 0.5.4 | CRAN (R 4.1.3) |
| ggstance | 0.3.6 | CRAN (R 4.1.3) |
| gld | 2.6.5 | CRAN (R 4.1.3) |
| globals | 0.16.1 | CRAN (R 4.1.3) |
| glue | 1.6.2 | CRAN (R 4.1.3) |
| googledrive | 2.0.0 | CRAN (R 4.1.2) |
| googlesheets4 | 1.0.1 | CRAN (R 4.1.3) |
| gower | 1.0.0 | CRAN (R 4.1.2) |
| gridExtra | 2.3 | CRAN (R 4.1.2) |
| gtable | 0.3.1 | CRAN (R 4.1.3) |
| hardhat | 1.2.0 | CRAN (R 4.1.3) |
| haven | 2.5.1 | CRAN (R 4.1.3) |
| highr | 0.9 | CRAN (R 4.1.3) |
| Hmisc | 4.7-1 | CRAN (R 4.1.3) |
| hms | 1.1.2 | CRAN (R 4.1.3) |
| htmlTable | 2.4.1 | CRAN (R 4.1.3) |
| htmltools | 0.5.3 | CRAN (R 4.1.3) |
| htmlwidgets | 1.5.4 | CRAN (R 4.1.2) |
| httpuv | 1.6.6 | CRAN (R 4.1.3) |
| httr | 1.4.4 | CRAN (R 4.1.3) |
| interp | 1.1-3 | CRAN (R 4.1.3) |
| ipred | 0.9-13 | CRAN (R 4.1.3) |
| jpeg | 0.1-9 | CRAN (R 4.1.1) |
| jsonlite | 1.8.2 | CRAN (R 4.1.3) |
| kableExtra | 1.3.4 | CRAN (R 4.1.3) |
| knitr | 1.40 | CRAN (R 4.1.3) |
| labeling | 0.4.2 | CRAN (R 4.1.1) |
| labelled | 2.10.0 | CRAN (R 4.1.3) |
| later | 1.3.0 | CRAN (R 4.1.2) |
| lattice | 0.20-45 | CRAN (R 4.1.1) |
| latticeExtra | 0.6-30 | CRAN (R 4.1.3) |
| lava | 1.6.10 | CRAN (R 4.1.2) |
| lifecycle | 1.0.3 | CRAN (R 4.1.3) |
| listenv | 0.8.0 | CRAN (R 4.1.2) |
| lme4 | 1.1-30 | CRAN (R 4.1.3) |
| lmom | 3.2 | CRAN (R 4.1.2) |
| lmtest | 0.9-40 | CRAN (R 4.1.3) |
| lubridate | 1.8.0 | CRAN (R 4.1.2) |
| magrittr | 2.0.3 | CRAN (R 4.1.3) |
| MASS | 7.3-58.1 | CRAN (R 4.1.3) |
| mathjaxr | 1.6-0 | CRAN (R 4.1.3) |
| Matrix | 1.5-1 | CRAN (R 4.1.3) |
| MatrixModels | 0.5-1 | CRAN (R 4.1.3) |
| memoise | 2.0.1 | CRAN (R 4.1.2) |
| MESS | 0.5.12 | CRAN (R 4.1.2) |
| meta | 7.0-0 | CRAN (R 4.1.2) |
| metadat | 1.2-0 | CRAN (R 4.1.2) |
| metafor | 4.2-0 | CRAN (R 4.1.2) |
| mice | 3.14.0 | CRAN (R 4.1.2) |
| mime | 0.12 | CRAN (R 4.1.1) |
| miniUI | 0.1.1.1 | CRAN (R 4.1.3) |
| minqa | 1.2.4 | CRAN (R 4.1.2) |
| mitools | 2.4 | CRAN (R 4.1.2) |
| modelr | 0.1.9 | CRAN (R 4.1.3) |
| mosaicCore | 0.9.2.1 | CRAN (R 4.1.3) |
| muhaz | 1.2.6.4 | CRAN (R 4.1.2) |
| multcomp | 1.4-20 | CRAN (R 4.1.3) |
| MuMIn | 1.46.0 | CRAN (R 4.1.2) |
| munsell | 0.5.0 | CRAN (R 4.1.2) |
| mvtnorm | 1.1-3 | CRAN (R 4.1.1) |
| nlme | 3.1-159 | CRAN (R 4.1.3) |
| nloptr | 2.0.3 | CRAN (R 4.1.3) |
| nnet | 7.3-18 | CRAN (R 4.1.3) |
| numDeriv | 2016.8-1.1 | CRAN (R 4.1.1) |
| pacman | 0.5.1 | CRAN (R 4.1.3) |
| parallelly | 1.32.1 | CRAN (R 4.1.3) |
| pillar | 1.8.1 | CRAN (R 4.1.3) |
| pkgbuild | 1.3.1 | CRAN (R 4.1.2) |
| pkgconfig | 2.0.3 | CRAN (R 4.1.2) |
| pkgload | 1.3.0 | CRAN (R 4.1.3) |
| png | 0.1-7 | CRAN (R 4.1.1) |
| polspline | 1.1.20 | CRAN (R 4.1.3) |
| polyclip | 1.10-4 | CRAN (R 4.1.3) |
| prettyunits | 1.1.1 | CRAN (R 4.1.2) |
| processx | 3.7.0 | CRAN (R 4.1.3) |
| prodlim | 2019.11.13 | CRAN (R 4.1.2) |
| profvis | 0.3.7 | CRAN (R 4.1.3) |
| progressr | 0.11.0 | CRAN (R 4.1.3) |
| promises | 1.2.0.1 | CRAN (R 4.1.2) |
| proxy | 0.4-27 | CRAN (R 4.1.3) |
| ps | 1.7.1 | CRAN (R 4.1.3) |
| purrr | 0.3.5 | CRAN (R 4.1.3) |
| quantreg | 5.94 | CRAN (R 4.1.3) |
| R6 | 2.5.1 | CRAN (R 4.1.3) |
| ragg | 1.2.3 | CRAN (R 4.1.3) |
| RColorBrewer | 1.1-3 | CRAN (R 4.1.3) |
| Rcpp | 1.0.9 | CRAN (R 4.1.3) |
| readr | 2.1.3 | CRAN (R 4.1.3) |
| readxl | 1.4.1 | CRAN (R 4.1.3) |
| recipes | 1.0.1 | CRAN (R 4.1.3) |
| remotes | 2.4.2 | CRAN (R 4.1.3) |
| renv | 1.0.1 | CRAN (R 4.1.2) |
| reprex | 2.0.2 | CRAN (R 4.1.3) |
| rlang | 1.0.6 | CRAN (R 4.1.3) |
| rmarkdown | 2.16 | CRAN (R 4.1.3) |
| rms | 6.3-0 | CRAN (R 4.1.3) |
| rootSolve | 1.8.2.4 | CRAN (R 4.1.2) |
| rpart | 4.1.16 | CRAN (R 4.1.3) |
| rstudioapi | 0.14 | CRAN (R 4.1.3) |
| Rttf2pt1 | 1.3.11 | CRAN (R 4.1.3) |
| rvest | 1.0.3 | CRAN (R 4.1.3) |
| sandwich | 3.0-2 | CRAN (R 4.1.3) |
| scales | 1.2.1 | CRAN (R 4.1.3) |
| sessioninfo | 1.2.2 | CRAN (R 4.1.2) |
| shiny | 1.7.2 | CRAN (R 4.1.3) |
| showtext | 0.9-5 | CRAN (R 4.1.3) |
| showtextdb | 3.0 | CRAN (R 4.1.3) |
| SparseM | 1.81 | CRAN (R 4.1.1) |
| stringi | 1.7.6 | CRAN (R 4.1.2) |
| stringr | 1.4.1 | CRAN (R 4.1.3) |
| survey | 4.1-1 | CRAN (R 4.1.2) |
| survival | 3.4-0 | CRAN (R 4.1.3) |
| svglite | 2.1.0 | CRAN (R 4.1.2) |
| sysfonts | 0.8.8 | CRAN (R 4.1.3) |
| systemfonts | 1.0.4 | CRAN (R 4.1.2) |
| tableone | 0.13.2 | CRAN (R 4.1.3) |
| textshaping | 0.3.6 | CRAN (R 4.1.3) |
| TH.data | 1.1-1 | CRAN (R 4.1.3) |
| tibble | 3.1.8 | CRAN (R 4.1.3) |
| tidylog | 1.0.2 | CRAN (R 4.1.3) |
| tidyr | 1.2.1 | CRAN (R 4.1.3) |
| tidyselect | 1.2.0 | CRAN (R 4.1.3) |
| tidyverse | 1.3.2 | CRAN (R 4.1.3) |
| timeDate | 4021.106 | CRAN (R 4.1.3) |
| tweenr | 2.0.2 | CRAN (R 4.1.3) |
| tzdb | 0.3.0 | CRAN (R 4.1.3) |
| urlchecker | 1.0.1 | CRAN (R 4.1.3) |
| usethis | 2.1.6 | CRAN (R 4.1.3) |
| utf8 | 1.2.2 | CRAN (R 4.1.2) |
| vcd | 1.4-10 | CRAN (R 4.1.3) |
| vctrs | 0.5.2 | CRAN (R 4.1.3) |
| viridisLite | 0.4.1 | CRAN (R 4.1.3) |
| visNetwork | 2.1.2 | CRAN (R 4.1.3) |
| webshot | 0.5.4 | CRAN (R 4.1.3) |
| withr | 2.5.0 | CRAN (R 4.1.2) |
| xfun | 0.33 | CRAN (R 4.1.3) |
| xml2 | 1.3.3 | CRAN (R 4.1.2) |
| xtable | 1.8-4 | CRAN (R 4.1.2) |
| yaml | 2.3.6 | CRAN (R 4.1.3) |
| zoo | 1.8-11 | CRAN (R 4.1.3) |
Code
reticulate::py_list_packages()%>%
knitr::kable(caption = "Python packages", format = "html",
col.names = c("Package", "Version", "Requirement"),
row.names = FALSE,
align = c("c", "l", "r", "r"))%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::scroll_box(width = "100%", height = "375px") | Package | Version | Requirement |
|---|---|---|
| absl-py | 2.1.0 | absl-py==2.1.0 |
| asttokens | 2.4.1 | asttokens==2.4.1 |
| astunparse | 1.6.3 | astunparse==1.6.3 |
| audioconverter | 2.0.3 | audioconverter==2.0.3 |
| autograd | 1.6.2 | autograd==1.6.2 |
| autograd-gamma | 0.5.0 | autograd-gamma==0.5.0 |
| beautifulsoup4 | 4.12.3 | beautifulsoup4==4.12.3 |
| Brotli | 1.1.0 | Brotli==1.1.0 |
| certifi | 2023.11.17 | certifi==2023.11.17 |
| cffi | 1.16.0 | cffi==1.16.0 |
| charset-normalizer | 3.3.2 | charset-normalizer==3.3.2 |
| clarabel | 0.9.0 | clarabel==0.9.0 |
| click | 8.1.7 | click==8.1.7 |
| cloudpickle | 3.0.0 | cloudpickle==3.0.0 |
| colorama | 0.4.6 | colorama==0.4.6 |
| comm | 0.2.1 | comm==0.2.1 |
| contourpy | 1.2.0 | contourpy==1.2.0 |
| cvxopt | 1.3.2 | cvxopt==1.3.2 |
| cvxpy | 1.5.2 | cvxpy==1.5.2 |
| cycler | 0.12.1 | cycler==0.12.1 |
| debugpy | 1.8.0 | debugpy==1.8.0 |
| decorator | 4.4.2 | decorator==4.4.2 |
| delete-chrome-history-py | 0.1.8 | delete-chrome-history-py==0.1.8 |
| easyocr | 1.7.1 | easyocr==1.7.1 |
| ecos | 2.0.13 | ecos==2.0.13 |
| editdistance | 0.8.1 | editdistance==0.8.1 |
| efficientnet | 1.0.0 | efficientnet==1.0.0 |
| essential-generators | 1.0 | essential-generators==1.0 |
| et-xmlfile | 1.1.0 | et-xmlfile==1.1.0 |
| executing | 2.0.1 | executing==2.0.1 |
| fancyimpute | 0.7.0 | fancyimpute==0.7.0 |
| ffmpeg | 1.4 | ffmpeg==1.4 |
| ffmpeg-python | 0.2.0 | ffmpeg-python==0.2.0 |
| filedir | 0.0.3 | filedir==0.0.3 |
| filelock | 3.13.1 | filelock==3.13.1 |
| flatbuffers | 24.3.25 | flatbuffers==24.3.25 |
| fonttools | 4.47.2 | fonttools==4.47.2 |
| formulaic | 1.0.1 | formulaic==1.0.1 |
| fsspec | 2023.12.2 | fsspec==2023.12.2 |
| future | 0.18.3 | future==0.18.3 |
| gast | 0.6.0 | gast==0.6.0 |
| git-filter-repo | 2.45.0 | git-filter-repo==2.45.0 |
| google-pasta | 0.2.0 | google-pasta==0.2.0 |
| graphviz | 0.20.3 | graphviz==0.20.3 |
| grpcio | 1.65.4 | grpcio==1.65.4 |
| gTTS | 2.5.1 | gTTS==2.5.1 |
| h5py | 3.11.0 | h5py==3.11.0 |
| idna | 3.6 | idna==3.6 |
| imageio | 2.34.2 | imageio==2.34.2 |
| imageio-ffmpeg | 0.5.1 | imageio-ffmpeg==0.5.1 |
| imgaug | 0.4.0 | imgaug==0.4.0 |
| iniconfig | 2.0.0 | iniconfig==2.0.0 |
| interface-meta | 1.3.0 | interface-meta==1.3.0 |
| ipykernel | 6.29.5 | ipykernel==6.29.5 |
| ipython | 8.20.0 | ipython==8.20.0 |
| jedi | 0.19.1 | jedi==0.19.1 |
| Jinja2 | 3.1.3 | Jinja2==3.1.3 |
| joblib | 1.4.0 | joblib==1.4.0 |
| jupyter_client | 8.6.0 | jupyter_client==8.6.0 |
| jupyter_core | 5.7.1 | jupyter_core==5.7.1 |
| keras | 3.4.1 | keras==3.4.1 |
| Keras-Applications | 1.0.8 | Keras-Applications==1.0.8 |
| keras-ocr | 0.9.3 | keras-ocr==0.9.3 |
| kiwisolver | 1.4.5 | kiwisolver==1.4.5 |
| knnimpute | 0.1.0 | knnimpute==0.1.0 |
| lazy_loader | 0.4 | lazy_loader==0.4 |
| libclang | 18.1.1 | libclang==18.1.1 |
| lifelines | 0.28.0 | lifelines==0.28.0 |
| llvmlite | 0.41.1 | llvmlite==0.41.1 |
| Markdown | 3.6 | Markdown==3.6 |
| markdown-it-py | 3.0.0 | markdown-it-py==3.0.0 |
| MarkupSafe | 2.1.4 | MarkupSafe==2.1.4 |
| matplotlib | 3.8.2 | matplotlib==3.8.2 |
| matplotlib-inline | 0.1.6 | matplotlib-inline==0.1.6 |
| mdurl | 0.1.2 | mdurl==0.1.2 |
| mido | 1.3.3 | mido==1.3.3 |
| ml-dtypes | 0.4.0 | ml-dtypes==0.4.0 |
| more-itertools | 10.2.0 | more-itertools==10.2.0 |
| moviepy | 1.0.3 | moviepy==1.0.3 |
| mpmath | 1.3.0 | mpmath==1.3.0 |
| multipledispatch | 1.0.0 | multipledispatch==1.0.0 |
| mutagen | 1.47.0 | mutagen==1.47.0 |
| namex | 0.0.8 | namex==0.0.8 |
| natsort | 8.4.0 | natsort==8.4.0 |
| nest-asyncio | 1.5.9 | nest-asyncio==1.5.9 |
| networkx | 3.2.1 | networkx==3.2.1 |
| ninja | 1.11.1.1 | ninja==1.11.1.1 |
| nose | 1.3.7 | nose==1.3.7 |
| numba | 0.58.1 | numba==0.58.1 |
| numexpr | 2.10.0 | numexpr==2.10.0 |
| numpy | 1.26.3 | numpy==1.26.3 |
| openai-whisper | 20231117 | openai-whisper==20231117 |
| opencv-python | 4.10.0.84 | opencv-python==4.10.0.84 |
| opencv-python-headless | 4.10.0.84 | opencv-python-headless==4.10.0.84 |
| openpyxl | 3.1.4 | openpyxl==3.1.4 |
| opt-einsum | 3.3.0 | opt-einsum==3.3.0 |
| optree | 0.12.1 | optree==0.12.1 |
| osqp | 0.6.5 | osqp==0.6.5 |
| packaging | 23.2 | packaging==23.2 |
| pandas | 2.2.0 | pandas==2.2.0 |
| pandas-flavor | 0.6.0 | pandas-flavor==0.6.0 |
| parso | 0.8.3 | parso==0.8.3 |
| patsy | 0.5.6 | patsy==0.5.6 |
| pillow | 10.2.0 | pillow==10.2.0 |
| platformdirs | 4.1.0 | platformdirs==4.1.0 |
| pluggy | 1.5.0 | pluggy==1.5.0 |
| polars | 1.9.0 | polars==1.9.0 |
| proglog | 0.1.10 | proglog==0.1.10 |
| prompt-toolkit | 3.0.43 | prompt-toolkit==3.0.43 |
| protobuf | 4.25.4 | protobuf==4.25.4 |
| psutil | 5.9.8 | psutil==5.9.8 |
| pure-eval | 0.2.2 | pure-eval==0.2.2 |
| pyarrow | 15.0.0 | pyarrow==15.0.0 |
| pyclipper | 1.3.0.post5 | pyclipper==1.3.0.post5 |
| pycparser | 2.22 | pycparser==2.22 |
| pycryptodomex | 3.20.0 | pycryptodomex==3.20.0 |
| pydotplus | 2.0.2 | pydotplus==2.0.2 |
| pydub | 0.24.1 | pydub==0.24.1 |
| Pygments | 2.17.2 | Pygments==2.17.2 |
| pyjanitor | 0.26.0 | pyjanitor==0.26.0 |
| PyMuPDF | 1.24.9 | PyMuPDF==1.24.9 |
| PyMuPDFb | 1.24.9 | PyMuPDFb==1.24.9 |
| pyparsing | 3.1.1 | pyparsing==3.1.1 |
| PyPDF2 | 3.0.1 | PyPDF2==3.0.1 |
| pyreadr | 0.5.0 | pyreadr==0.5.0 |
| pytesseract | 0.3.10 | pytesseract==0.3.10 |
| pytest | 8.3.1 | pytest==8.3.1 |
| python-bidi | 0.6.0 | python-bidi==0.6.0 |
| python-dateutil | 2.8.2 | python-dateutil==2.8.2 |
| pytube | 15.0.0 | pytube==15.0.0 |
| pytube3 | 9.6.4 | pytube3==9.6.4 |
| pytz | 2023.3.post1 | pytz==2023.3.post1 |
| pywin32 | 306 | pywin32==306 |
| PyYAML | 6.0.1 | PyYAML==6.0.1 |
| pyzmq | 25.1.2 | pyzmq==25.1.2 |
| qdldl | 0.1.7.post1 | qdldl==0.1.7.post1 |
| regex | 2023.12.25 | regex==2023.12.25 |
| requests | 2.32.3 | requests==2.32.3 |
| rich | 13.7.1 | rich==13.7.1 |
| rpy2 | 3.5.16 | rpy2==3.5.16 |
| scikit-image | 0.24.0 | scikit-image==0.24.0 |
| scikit-learn | 1.3.2 | scikit-learn==1.3.2 |
| scikit-survival | 0.22.2 | scikit-survival==0.22.2 |
| scipy | 1.11.4 | scipy==1.11.4 |
| scs | 3.2.6 | scs==3.2.6 |
| seaborn | 0.13.2 | seaborn==0.13.2 |
| semantic-version | 2.10.0 | semantic-version==2.10.0 |
| setuptools-rust | 1.8.1 | setuptools-rust==1.8.1 |
| shapely | 2.0.5 | shapely==2.0.5 |
| six | 1.16.0 | six==1.16.0 |
| soupsieve | 2.5 | soupsieve==2.5 |
| SpeechRecognition | 3.10.1 | SpeechRecognition==3.10.1 |
| spyder-kernels | 2.5.2 | spyder-kernels==2.5.2 |
| stack-data | 0.6.3 | stack-data==0.6.3 |
| statsmodels | 0.14.1 | statsmodels==0.14.1 |
| sympy | 1.12 | sympy==1.12 |
| target | 0.0.11 | target==0.0.11 |
| tensorboard | 2.17.0 | tensorboard==2.17.0 |
| tensorboard-data-server | 0.7.2 | tensorboard-data-server==0.7.2 |
| tensorflow | 2.17.0 | tensorflow==2.17.0 |
| tensorflow-intel | 2.17.0 | tensorflow-intel==2.17.0 |
| tensorflow-io-gcs-filesystem | 0.31.0 | tensorflow-io-gcs-filesystem==0.31.0 |
| termcolor | 2.4.0 | termcolor==2.4.0 |
| threadpoolctl | 3.4.0 | threadpoolctl==3.4.0 |
| tifffile | 2024.7.24 | tifffile==2024.7.24 |
| tiktoken | 0.5.2 | tiktoken==0.5.2 |
| torch | 2.4.0 | torch==2.4.0 |
| torchaudio | 2.4.0 | torchaudio==2.4.0 |
| torchvision | 0.19.0 | torchvision==0.19.0 |
| tornado | 6.4 | tornado==6.4 |
| tqdm | 4.66.1 | tqdm==4.66.1 |
| traitlets | 5.14.1 | traitlets==5.14.1 |
| translator | 0.0.9 | translator==0.0.9 |
| typing_extensions | 4.9.0 | typing_extensions==4.9.0 |
| tzdata | 2023.4 | tzdata==2023.4 |
| tzlocal | 5.2 | tzlocal==5.2 |
| urllib3 | 2.1.0 | urllib3==2.1.0 |
| validators | 0.33.0 | validators==0.33.0 |
| watchdog | 3.0.0 | watchdog==3.0.0 |
| wcwidth | 0.2.13 | wcwidth==0.2.13 |
| websockets | 12.0 | websockets==12.0 |
| Werkzeug | 3.0.3 | Werkzeug==3.0.3 |
| whisper | 1.1.10 | whisper==1.1.10 |
| wrapt | 1.16.0 | wrapt==1.16.0 |
| xarray | 2024.1.1 | xarray==2024.1.1 |
| youtube-dl | 2021.12.17 | youtube-dl==2021.12.17 |
| yt-dlp | 2024.7.9 | yt-dlp==2024.7.9 |
Output
Code
invisible("Adicionales en E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/para_pres_cum_mean.R")
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
"E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
"G:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")
save.image(paste0(folder_path,"an_grant_23_24_4.RData"))
invisible("ültima vs. estable")
#load(paste0(folder_path,"1_an_grant_23_24_4.RData"))